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JET FORMATION IN BLACK HOLE ACCRETION SYSTEMS II: NUMERICAL MODELS 

Jonathan C. McKinney 1 

June 16, 2005 

ABSTRACT 

In a companion theory paper, we presented a unified model of jet formation. We suggested that primarily two 
types of relativistic jets form near accreting black holes: a potentially ultrarelativistic Poynting-dominated jet 
and a Poynting-baryon jet. We showed that, for the collapsar model, the neutrino-driven enthalpy flux (classic 
fireball model) is probably dominated by the Blandford-Znajek energy flux, which predicts a jet Lorentz factor 
of r ~ 100- 1000. We showed that radiatively inefficient AGN, such as M87, are synchrotron-cooling limited 
to r ~ 2- 10. Radiatively efficient x-ray binaries, such as GRS1915+105, are Compton-drag limited to Y < 2, 
but the jet may be destroyed by Compton drag. However, the Poynting-baryon jet is a collimated outflow with 
r ~ 1-3. Here we present general relativistic hydromagnetic simulations of black hole accretion with pair 
creation used to simulate jet formation in GRBs, AGN, and x-ray binaries. Our collapsar model shows the 
development of a patchy "magnetic fireball" with typically Y ~ 100- 1000 and a Gaussian structure. Temporal 
variability of the jet is dominated by toroidal field instabilities for > 10 2 gravitational radii. A broader Poynting- 
baryon jet with r ~ 1.5 could contribute to a supernova. 

Subject headings: accretion disks, black hole physics, galaxies: jets, gamma rays: bursts, X-rays : bursts, 
supernovae: general, neutrinos 



1. INTRODUCTION 

General relativistic magnetohydrodynamics (GRMHD) is 
the black hole mass-invariant (nonradiative) physics com- 
monly used to describe black hole accretion systems. Such 
systems often exhibit jets. However, the observed jet prop- 
erties, such as co llimation and spee d, are not uniform be- 
tween systems. In McKinnev (2005b), our goal was to deter- 
mine the unifying, or minimum number of, pieces of physics 
that would explain jet formation associated with gamma-ray 
bursts (GRBs), active galactic nuclei (AGN), and x-ray bi- 
naries. Similar ideas have been explored by other authors 
(Ghise llini & Celottil l200l iGhisellinil 120031 iMelerl [2003). 
The goal of McKinnev (2005b) was to understand jet forma- 
tion and explain the origin of the energy, composition, colli- 
mation, and Lorentz factor. 

This paper impl ements some of the theoretical ideas of 
IMcKinnevi ( 2005b ) into a GRMHD numerical model. In par- 
ticular, we study the collapsar model with pair creation and 
an effective model for Fick diffusion of neutrons that con- 
taminate the jet with an electron-proton plasma. The generic 
model parameters also allow the numerical model to be inter- 
preted as a model for AGN systems such as M87. We also 
comment on the applicability of the model to x-ray binary 
systems. 

§ [2] summarizes the proposed unified model to explain jet 
formation in all blac k hole accretion sys tems, where more de- 
tails are provided in McKinnev (2005b). 

§[3]presents the results of a fiducial numerical model corre- 
sponding to the collapsar model and describes the jet structure 
and formation. The relevance of the presented model to AGN 
and x-ray binaries is discussed. The jet is found to have some 
piece-wise self-similar features, and curve fits are given for 
use by other modellers. Acceleration of the GRB jet is found 
to occur over a large range in radii rather than occurring close 
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to the bl ack hole. T he introductory estimates of the Lorentz 
factor in McKinnev ( 2005bj are refined based upon these nu- 
merical models. The jet characteristic structure, such as the 
formation of a double transonic (transfast) flow, is discussed. 

§ 0] discusses the results and their possible implications. 
The results are compared to similar investigations, and the 
limitations of the models presented are discussed. 

§[5]summarizes the key points. 

In appendix |A| we give the GRMHD equations of motion 
solved when accounting for pair creation. See also lMcKinnevI 
(2005b) for relevant discussions. These are the equations nu- 
merically solved that give the results discussed in section [3] 
These are straightforward and how one handles the details 
of the pair creation ends up not making any qualitative (or 
much quantitative) difference in the results, hence why the 
material is in the appendix. That is, the energetics of the 
Poynting-dominated jet are dominated by the electromagnetic 
field. AppendixlElsumm arizes the well-known characteristics 
of the GRMHD equations and other surfaces of physical in- 
terest used in the discussion in section|3] 

2. GRMHD PAIR INJECTION MODEL OF JET FORMATION 

In IMcKinnevi i2005b). we showed that one can identify a 
small subset of physics that can explain the jet energy, com- 
position, collimation, and Lorentz factor for all black hole ac- 
cretion systems. 

One key idea of that paper is that the terminal Lorentz fac- 
tor is determined by the toroidal magnetic energy per unit pair 
mass density energy near the location where pairs can escape 
to infinity (beyond the so-called "stagnation surface"). For 
GRBs, neutron diffusion is crucial to explain (and limit) the 
Lorentz factor and explain the baryon-pollution or baryon- 
contamination problem. For AGN and x-ray binaries, since 
a negligible number of baryons cross the field lines, pair- 
loading from radiative annihilation is crucial to determine the 
Lorentz factor of the Poynting-dominated jet since this deter- 
mines the rest-mass flux and density in the jet. 
E Figure ^ shows the basic picture for GRB systems, while 
figure [2] shows the basic picture for AGN and x-ray binary 
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FIG. 1. — Schematic of pair-production model and subsequent magnetic 
fireball formation for GRB disks. Fireball is extremely optically thick. Below 
a stagnation surface, pairs are accreted by the black hole and so do not load 
the jet. Here Rg = GM/c 2 . 
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FIG. 2. — Schematic of pair-production model and subsequent shock- 
heating and emission. AGN jet is optically thin and emits nonthermal and 
thermal synchrotron, while x-ray binary jet can be marginally optically thick 
and emit via self-absorbed synchrotron and by severe Compton drag. Severe 
Compton drag can lead to destruction of the Poynting-dominated jet. 



systems. An accreting, spinning black hole creates a mag- 
netically dominated funnel region around the polar axis. The 
rotating black hole drives a Poynting flux into the funnel re- 
gion, where the Poynting flux is associated with the coiling of 
poloidal magnetic field lines into toroidal magnetic field lines. 
The ideal MHD approximation holds very well, and so only 
neutral particles such as photons, neutrinos, and neutrons can 
cross the field lines and load the Poynting-dominated jet that 
emerges. 

The accretion disk emits neutrinos in a GRB model (7-ray 
and many soft photons for AGN and x-ray binaries) that an- 
nihilate and pair-load the funnel region within some "injec- 
tion region." For GRB systems, neutrons Fick-diffuse across 
the field lines and collisionally decay into an electron-proton 
plasma. 

Many pairs (any type) are swallowed by the black hole, but 
some escape if beyond some "stagnation surface," where the 
time-averaged poloidal velocity is zero and positive beyond. 



Pairs beyond the stagnation surface are then accelerated by 
the Poynting flux in a self-consistently generated collimated 
outflow. In the electromagnetic (EM) jet, the acceleration pro- 
cess corresponds to a gradual uncoiling of the magnetic field 
and a release of the stored magnetic energy that originated 
from the spin energy of the black hole. 

One key result of this paper is that the release of magnetic 
energy need not be gradual once the toroidal field dominates 
the poloidal field, in which case pinch (and perhaps kink) in- 
stabilities can occur and lead to a nonlinear coupling (e.g. a 
shock) that converts P oynting flux into enthalpy flux ( Eichlei 
ll993t lBegerman fl998l) . In the proposed GRB model, this con- 
version reaches equipartition and the jet becomes a "magnetic 
fireball," where the toroidal field instabilities drive large vari- 
atio ns in the jet Lorentz factor and jet luminosity. 

In McKinnev (2005b), we showed that in AGN systems, 
nonthermal synchrotron from shock-accelerated electrons and 
some thermal synchrotron emission releases the shock en- 
ergy until the synchrotron cooling times are longer than the 
jet propagation time. For AGN, jet acceleration is negligi- 
ble beyond the extended shoc k zone, as sugg ested for blazars 
beyond the "blazar zone" ( Sikora et al. 2005). In x-ray binary 
systems, the shock is not as hot and also unlike in the AGN (at 
least those like M87) case the jet can be optically thick. Thus 
these x-ray binary systems self-absorbed synchrotron emit if 
they survive Compton drag. 

For all these systems, at large radii patches of energy flux 
and variations in the Lorentz factor develop due to toroidal 
instabilities. These patches in the jet could drive internal 
shocks and at large radii they drive external shocks with the 
surrounding medium. The EM jet is also surrounded by a 
mildly relativistic matter coronal outflow/jet/wind, which is 
a material extension of the corona surrounding the disk. This 
Poynting-baryon, coronal outflow collimates the outer edge of 
the Poynting-dominated jet, which otherwise internally col- 
limates by hoop stresses. The luminosity of the Poynting- 
baryon jet is determined, like the Poynting-dominated jet, by 
the mass accretion rate, disk thickness, and black hole spin. 

In this paper, this unified model is studied numerically 
using axisymmetric, nonradiative, GRMHD simulations to 
study the self-consistent process of jet formation from black 
ho le accretion systems. These s imulations extend the work 
of McKinnev & Gammie ( 2004) by including pair creation 
(and an effective neutron diffusion for GRB-type systems) to 
self-consistently treat the creation of jet matter, investigating 
a larger dynamic range in radius, and presenting a more de- 
tailed analysis of the Poynting-dominated jet structure. 

2.1. Units and Notation 

The units in this paper have GM = c = 1, which sets the 
scale of length (r g = GM/c 2 ) and time (f g = GM/c 1 ). The 
mass scale is determined by setting the (model-dependent) 
observed (or inferred for GRB-type systems) mass accretion 
rate [Mo[g s~ 1 ]) equal to the accretion rate through the black 
hole horizon as measured in a simulation. So the mass is 
scaled by the mass accretion rate (Mo) at the horizon (r = 
r H = r L (\ + y/l-j 2 )), such that p iUisk = M [r = r H ]t g /r 3 g and 

the mass scale is then just m = pojiskr 3 = Mo[r = rn]t g . Un- 
less explicitly stated, the magnetic field strength is given in 
Heaviside-Lorentz units, where the Gaussian unit value is ob- 
tained by multiplying the Heaviside-Lorentz value by \/4n. 

The value of poji S k can be determined for different sys- 
tems. For example, a collapsar model with M = 0. IMqs" 1 and 
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gem . M87 has a mass 
yr" 1 and a black ho le mass 
Reynolds et aTlH996l) giving 



M m 3M Q , then p 0i d isk « 3.4 x 10 1 
accretion rate of M ~ 10~ 2 Mr. 
of M « 3 x 10 9 M Q (lHolll999t 

Po.rf/sjt ~ 10~ 16 gcm" 3 . GRS 1915+105 has a mass accretion 
rate of M ~ 7 x lO^M^yr" 1 ilMirabel & Rodriguez! 11994 
iMirabel & Rodriguez! fl 9991: Fender & Bellonil 120041) with a 
mass o f M ~ 14Mq l lGreiner et al.120011) . but see lKaiser et al.l 
(2004). This gives /Oo,rfM ~ 3 x 10~ 4 gcm~ 3 . This disk density 
scales many of the results of the paper. The disk height (H) to 
radius (R) ratio is written as H/R. 

The GRMHD notation follows MTW. For example, the 
4-velocity components are m p (contravariant) or (covari- 
ant). For a black hole with angular momentum J = jGM 2 /c, 
j = a/Mis the dimensionless Kerr parameter with -1 <j< 1. 
The rest-mass density is given as po, internal energy density 

as u, magnetic field 3-vector as B' = F , where F is the 
dual of the Faraday tensor. The contravariant metric com- 
ponents are g^ u and covariant components are g^, where 
beyond the frame-dragging of the black hole the metric in 
Boyer-Lindquist coordinates is approximately diagonal such 
that an orthonormal contravariant component is ifi = jg^/ 
for any spatial component of m m . The comoving energy den- 
sity is b 2 /2, where £> M is the comoving magnetic field. See 
Mc Kinnev & Gammid ( 12004 for details. 

The Lorentz factor of the jet can be measured either as the 
current time-dependent value, or, using information about the 
GRMHD system of equations, one can estimate the Lorentz 
factor at large radii from fluid quantities at small radii. The 
Lorentz factor as measured by a static observer at infinity is 

r = M '" = M 'v^; (i) 

in Boyer-Lindquist coordinates, where no static observers 
exist inside the ergosphere. This is as opposed to W = 
u'y/—l /g", which is the Lorentz factor as measured by the 
nor mal observer as used by most numerical relativists. 

In lMcKinnevl (l2005b) we show that the Lorentz factor and 
(^-velocity at large distances are 

Y 00 =E = -hu t + m F B 4> (2) 

ui=L = hu 4> + ^>B 4> , (3) 

where h = (po+u g +p)/ po is the specific enthalpy, $ is the con- 
served magnetic flux per unit rest-mass flux, dp is the con- 
served field rotation frequency, B^ is the covariant toroidal 

magnetic field, and = M0 jOO . Here E and L simply rep- 
resent the conserved energy and angular momentum flux per 
unit rest-mass flux. Notice the matter and electromagnetic 
pieces are separable, such that 



r _ r (MA) , r (EM) 



<p,oo <p,oo ' 



(4) 
(5) 



See lMcKinnevl (f2005b) for more details. 



3. NUMERICAL EXPERIMENTS 

This section presents the results of GRMHD numeri- 
cal models with pair creation to self-consistently describe 
the Poynting-dominated and Poynting-baryo n jet formation 
proces s. Our numerical scheme is HARM ( Gam mie et al.l 
2003a), a conservative, shock-capturing scheme for evolv- 
ing the equations of general relativistic MHD. Compared to 
the original HARM, the inversion of conserved quantities to 
primitive variables is performed by using a new faster and 



more robust two-dimensional non-linear solver ( Noble et alJ 
120051) . The new H ARM also uses a parabolic interpola- 
tion scheme dColellall 19841) rather than a linear interpolation 
scheme. The new HARM also u ses an opt imal TVD third 
order Runge-Kutta time stepping (IShul 11997) rather than the 
mid-point method. For the problems under consideration, the 
parabolic interpolation and third order time stepping method 
reduce the truncation error significantly, including magneti- 
cally dominated regions with b 2 /po^> 1 . 

Notice that no explicit reconnection model is included. 
However, HARM checks the effective resistivity by measur- 
ing the rest-mass flux across field lines. For the boundary be- 
tween the coronal wind and the Poynting-dominated jet, the 
total mass flux across the field line is negligible compared to 
the rest-mass flux along the field line and the pair creation 
rate. An unresolved model would load field lines with rest- 
mass that crosses field li nes and no t properly represent any 
physical resistivity ( Mc Kinne\l2004 . 

3.1. Computational Domain 

The computational domain is axisymmetric, with a grid that 
typically extends from r ln = 0.98r#, where r# is the horizon, 
to r our = I0 4 r g , and from 8 = to 8 = tt. Full 3D calcula- 
tions with this dynamic range are not possible with today's 
computers. Our numerical integrations are carried out on a 
uniform grid in so-called "modified KS" (MKS) coordinates: 
xo,x\ ,X2,X3, where Xq = t [KS], X3 = </>[KS]. The radial coordi- 
nate is chosen to be 



r = /?o + e x < 



(6) 



where Rq is chosen to concentrate the grid zones toward the 
event horizon (as Rq is increased from to r#) and n r is con- 
trols the enhancement of inner to outer radial regions. For 
studies where the disk and jet interaction is of primary inter- 
est, Ro = and n r = 1 are chosen. For studies where in addition 
the far-field jet is of interest, Ro = -3 and n r = 10 are chosen. 
The 8 coordinate is chosen to be 



: nx 2 + ^(l-h(r)) sin(27rx 2 ), 



(7) 



where h(r) is used to concentrate grid zones toward the equa- 
tor (as h is decreased from 1 to 0) or pole (as h is increased 
from 1 to 2). The jet at large radii is resolved together with 
the disk at small radii using 

h(r) = 2-Qj(r/rojr' g ' (8) 

with the parameters of Qj = 1.3 — 1.8, r§j = 2.8, rij = 0.3, 
r\j = 20, r2j = 80, and gj = gj(r) = 5 + ^atani^ 1 ) is used to 
control the transition between inner and outer radial regions. 
To perform convergence tests, Qj and the overall radial and 
8 resolution are varied. An alternative to this fixed refine- 
ment of the jet and disk is an adaptive refinement (see, e.g., 
IZhang & M acFadv enl2005T) . 

3.2. Initial Conditions and Problem Setup 

All the experiments evolve a weakly magnetized torus 
around a Kerr black hole in axisymmetry. The focus of our 
numerical investigation is to study a high resolution model of 
the collapsar model. Other simulations were performed, and 
the fiducial model's relevance to AGN and x-ray binary sys- 
tems is discussed at the end of the section. A generic model 
is used so the results mostly can be applied to any black hole 
system. A black hole spin of j = 0.9375 is chosen, but this 
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produces similar results to models with 0.5 < j < 0.97, which 
includes most black hole accretion systems. 

The initial conditions consist of an equilibrium 
torus which is a "donut" of plasma with a black 
hole at the center (iFishb one & Moncriell 119761: 
Abramowicz. Jaroszinski. & Sikora 1978). The donut is 
supported again st gravity by centrifugal and pressure forces. 
The solutions of Fishbone & Moncrief ( 1976) corresponding 
to u'u^ = const, are used. The initial inner edge of the torus 
is set at r C£ / gc = 6. The equation of state is a gamma-law 
ideal gas with 7 = 4/3, but other 7 lead to similar results 
(McKinnev & Gammie 2004). This donut is a solution to the 
axisymmetric stationary equations, but the donut is unstable 
to gl obal nonaxisymmetric modes ( Papaloizou & Pringle 
1983). However, when the donut is embedded in a weak 
magnetic field, the magnetorotational instability (MRI) 
dominates those hydrodynamic modes. Small perturbations 
are introduced in the velocity field to seed the instability. The 
models have u'u^ = 4.281, the pressure maximum is located 
at r max = 12, the inner edge at (r, 9) = (6,7r/2), and the outer 
edge at (r,(9) = (42, tt/2). 

This is consistent with the accretion disk th at is expected to 
form in the collapsar model (MacFadv en & WooslevlH999T) . 
A GRB in the collapsar model is formed after the collapse 
of a massive star as a black hole with Mbh ~ 3M Q and spin 
j ~ 0.75-0.94 accretes at Mo ~ 0.1M©/s from a torus of mat- 
ter within r < 200km that has a hei ght (H) to radius (R) ra- 
tio ("th ickness") of H/R ~ 0.2-0.4. MacFa dven & Wooslevl 
( 1999) find that an energetic outflow develops in the polar re- 
gion where p ~ 10 7 gcm" 3 near the horizon, p ~ 10 5 gcm" 3 
at r ~ 200 km, and p ~ 10 3 ~ 4 gcm~ 3 at r ~ 2000km (see 
model 14A in MacFadven & Wooslev 1999). In the par- 
lance of collapsar models, our u ! u^ model corresponds to a 
/us = //(10 16 cm 2 s~ 1 ) of /i f, ~ 5, within the range suggested 
bv iMacFadven & Wooslevl (fl999) to produce a GRB. Also, 
the H/R^ 0.2-0.4 as suggested th at forms according to neu- 
trino emission models (Kohri & Mineshige 2002; Koh ri et aTl 
12001 . 

For radiatively inefficient AGN and x-ray binaries this is 
simply a reasonable approximate model for the inner-radial 
accretion flow. While a slightly thicker disk (H/R^ 0.9) may 
be more appropriate for radiatively inefficient flows, this is 
not expected to significantly affect these results. A study of 
the effect of disk thickness, especially thin disks, is left for 
future work. 

The orbital period at the pressure maximum 27r(a + 
(fmax/ r g) 3 2 )tg — 267 t g , as measured by an observer at infin- 
ity. The model is run for At = 1.4 x 10 4 f g , which is about 
52 orbital periods at the pressure maximum and about 1150 
orbital periods at the black hole horizon. For the collapsar 
model this is only ~ 0.2 seconds, and at the spin chosen would 
correspond to the late phase of accretion when the Poynting 
flux reaches its largest magnitude. For the AGN M87 this cor- 
responds to ~ 7 years. For the x-ray binary GRS 1915+105 
this corresponds to ~ 1 second. 

Notice that since the model is axisymmetric, disk turbu- 
lence is not sustained after about t ~ 3000? ? ; that is, the anti- 
dynamo theorem prevails (Cowlina ll934l) . However, while 
this affects the disk accretion, this does not affect the evolu- 
tion of the Poynting-dominated jet. That is, from the time of 
turbulent accretion to "laminar" accretion, the funnel region 
is mostly unchanged. Indeed, the far-field jet that has already 
formed is causally disconnected from the region where the 



accretion disk would still be turbulent. 

For the collapsar model, at the quasi-steady state mass ac- 
cretion rate of 0. 1M Q /s the disk will last for t < 0.42s. Thus, 
this model requires a more extended disk or a fresh supply 
of plasma into the disk from the surrounding stellar enve- 
lope throug h an "accretion shock" to gen erate a long dura- 
tion GRBs (MacFad ven & Wooslevl fl999T) . However, since 
any such system modelled is in quasi-steady state, the results 
here do not strongly depend on the mass of the disk as long as 
the disk is not too massive, which would require taking into 
account the self-gravity of the disk. 

The numerical resolutio n of these models is 5 12 x 256 com- 
pared to 456 2 in Mc Kinnev & Gammid (120041) . However, due 
to the enhanced 9 grid, the resolution in the far-field jet region 
is ~ 10 times larger. Also, with the use of a parabolic interpo- 
lation scheme, the overall resolution is additionally enhanced. 
Compared to our previous model this gives us an effective 9 
resolution of ss 9000. 

As with lMcKinnev &~Ga mmie (2004), into the initial torus 
is put a purely poloidal magnetic field. The field can be de- 
scribed using a vector potential with a single nonzero com- 
ponent oc MAX(/9o//3o.mfl.v _ 0.2,0) The field is therefore 
restricted to regions with po/po,max > 0.2. The field is nor- 
malized so that the minimum ratio of gas to magnetic pres- 
sure is 100. The equilibrium is therefore only weakly per- 
turb ed by the magnetic field . It is, however, n o longer sta- 
ble llBalbus & Hawlevlll991l iGammid 120041) . iHirose et all 
( 120041) : iMcKinnev & Gammid (12004 show that no initial 
large-scale net vertical field is necessary, since a large-scale 
poloidal field is self-consistently generated. The field con- 
nects the black hole horizon (as observed by a physical ob- 
server) to large distances. Apart from the existence of an over- 
all poloidal sign, the results are insensitive to the details of the 
initial magnetic field geometry for any physically motivated 
geometry (McK innev & Gammiel20 04). 

3.3. Floor vs. Pair Creation Model 

Numerical models often "model" the injection physics in 
the Poynting-dominated jet by employing a so-called floor 
model. This model forces a minimum on the rest-mass den- 
sity {pfi) and internal energy density («//), which are usually 
set to several orders of magnitude lower than the disk density. 
Actually, however, floor models have always been treated as 
a purely numerical invention in order to avoid numerical ar- 
tifacts associated with the inability to numerically solve the 
equations of motion when the density is low. Indeed, the idea 
was one should convergence test by gradually lowering pji 
and Ufi. This is perhaps a reasonable for numerical study of 
stars, but accretion flows are so hot (baryons with T > \Q m K 
in the thick disk state near the black hole) that pairs are cre- 
ated when neutrinos annihilate (or when photons annihilate 
for x-ray binary and AGN systems). 

In previous numerical models of Poynting jets, the funnel 
region is always completely evacuated, and floor-models nec- 
essarily generate mass at least where the poloidal velocity 
u p = at the stagnation surface. That is, matter inside the stag- 
nation surface falls into the black hole, while matter outside 
it is ejected as part of the jet. Indeed, arbitrary floor-models 
violate the ideal-MHD condition far away from the black hole 
where no pairs should be produced. For example, if a numer- 
ical model uses pji ~ Const., then this leads to a completely 
unphysical model of pair creation and the ideal-MHD approx- 
imation is violated for the entire length of the jet in the funnel. 

While pair creation has so far been ignored by those doing 
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numerical models of jets from black hole accretion systems, it 
has been shown that the properties of the Poynting-dominated 
jet are almost completely determine d by the pair creation 
physics, which thus c annot be ignored (Phinney: 1983; Punslv 
199lt lLevinsonl2005l) . 

Our collapsar model accounts for pair creation and Fick dif- 
fusion of neutrons. In this paper, there is an injected rest- 
mass, enthalpy, and momentum. For the fiducial collapsar- 
like model, rest-mass, enthalpy, and momentum are injected 
with energy at infinity fractions of, res pectively, f n = 0.05 , 
//, = 0.45, and f m = 0.5 as described in McKinney (2005b). 
The total energy at infinity injected is determined by neutrino 
annihilation rates and Fick diffusion rates. The energy injec- 
tion rate from neutrino annihi l ation is based upon viscous disk 
models (Popha m et aTlll999t iMacFadven & Wooslevlll999l) . 
where the rate of injection is determined by choosing a vis- 
cosity coefficient a = 0.01 and mass accretion rate Mq = 
O.IM0S" 1 for the collapsar model. As discussed in that pa- 
per, the rest-mass is dominated by electron-positron pairs only 
very close to the black hole. 

The value of a and Mq determine the energy injected in 
proportion to the disk density per unit light crossing time (t g ). 
Thus any system with comparable normalized energy injec- 
tion rate would follow similar results. For example, the M87 
model is similar, where there the electron-positron pairs do 
not annihilate and represent the true density in the jet. In the 
M87 model there is no Fick diffusion, yet the dimensionless 
model parameters are approximately the same. 

The momentum direction is chosen by setting uf n j = i& and 

u fnj = 0' so mat fm determines u r m j. It turns out that the choice 
of /;,, f m , and the injected momentum direction very weakly 
determine the jet structure for r > 6r g . This is because a large 
fraction of the pairs are lost to the black hole and magnetic 
forces, rather than the injected momentum, dominate the mo- 
tion of the plasma near the black hole. 

Beyond this injection region, the rest-mass is dominated 
by electron-protons created in collision events from neutrons 
that Fick diffuse across the field lines into the jet region. For 
the collapsar model this coincidentally corresponds to the in- 
jected rest-mass injected in this model. Since the electron- 
positron pair annihilation only gives an additional ~ 10% of 
internal energy, then pair annihilation need not be considered, 
and the rest-mass injected well-models the rest-mass injected 
as an electron-proton plasma from Fick-diffused neutrons that 
collisionally decay. 

This physical model of the injected particles is augmented 
when the rest-mass or internal energy density reaches very 
low values. If the density drops below some "floor" value, 
then the density is returned to the floor value. The floor model 
chosen has pji = 10"' 1 r~ 2n pojisk and Ufi = 10~ 9 r~ 2 ' 7 pojisk- 
HARM keeps track of how often the density goes below the 
floors and how this modifies the conservation of mass, en- 
ergy, and angular momentum. The floor model contribution 
is negligible compared to the physical injection model at all 
spatial locations and at all times. The coefficient of the floor 
was chosen to be close to the resulting density near the black 
hole horizon. It was not chosen arbitrarily small in order to 
avoid numerical difficulties in integrating the equations when, 
rarely, the floor is activated in regions of large magnetic en- 
ergy density per unit rest-mass density. 

3.4. Boundary Conditions 



Two different models are chosen for the outer boundary 
condition. One model is called an "outflow" boundary con- 
dition, for which all primitive variables are projected into the 
ghost zones while forbidding inflow. The other model is to 
inject matter at some specified rate at the outer boundary, un- 
less there is outflow. This would correspond to a Bondi-like 
infall for AGN and x-ray binaries, or would correspond to 
the collapsing envelope of a massive star for collapsar mod- 
els. In particular, unless there is outflow, the outer bound- 
ary is set to inject mass at the free-fall rate, with a density 
of pa = I0~ 4 r~ 3 ' 2 polish and u = 10" 6 r" 5 ' 2 po^isk and no angular 
momentum. Presupernova models suggest that the infalling 
matter is at about 30% the free -fall rate we have chosen above 
( MacF adven & Wooslevll 99*91) . but this difference is unlikely 
to significantly affect the jet formation. Also, for their collap- 
sar model, the density structure in the equatorial plane varies 
between po oc r~ 3 / 2 to po oc r~ 2 and the internal energy is 
u oc r~ 5//2 to u oc r~ 2 7 (MacFadven & Woosle vll999l) . This is 
similar to our simplified model, which happens to also model 
a Bondi accretion for AGN and x-ray binaries. Thus the re- 
sults are indicative of GRBs, AGN, and x-ray binaries. The 
outer grid radius corresponds to about 20 presupernova core 
radii dWooslev & Weaver! 1995) or ~ 10 10 km or about 1 /10th 
the entire star's radius. 

3.5. Fiducial Model 

The overall character of the accret ion flow is u nchanged 
compa red to the descriptions given in McKi nnev & Ga mmie 
(2004). The disk enters a long, quasi-steady phase in which 
the accretion rates of rest-mass, angular momentum, and en- 
ergy onto the blac k hole fluctuate around a well-d efined mean. 
Meanwhile, as in McKinnev & Gammie (2004), a Poynting- 
dominated jet and Poynting-baryonjet (coronal outflow) have 
formed. The Poynting-dominated jet forms once the ram pres- 
sure of the funnel material is lower than the toroidal mag- 
netic pressure. Afterwards the baryons are spread apart in the 
launch of a magnetic tube filled with pairs. This occurs within 
t < 500f g . 

The Poynting-dominated jet forms as the differential rota- 
tion of the disk and the frame-dragging of the black hole in- 
duce a significant toroidal field that launches materia l away 
from the black hole by the same force described in McKinnev 
(20053. 

A coronal outflow is also generated between the disk and 
Poynting-dominated jet. In this model the coronal outflow 
has Too ~ 1.5. The coronal-funnel boundary contains shocks 
with a sonic Mach number of M s ~ 100. The inner-radial 
interface between the disk and corona is a site of vigorous 
reconnection due to the magnetic buoyancy and convective 
instabilities present there. These two parts of the corona are 
about 100 times hotter than the bulk of the disk. Thus these 
coronal components are a likely sites for Comptonization and 
nonthermal particle acceleration. 

Figure [3] and figure 0] show the final time (f = 1.4 x 
lO 4 ^) log of rest-mass density and magnetic field projected 
on the Cartesian z vs. x plane. For the purposes of 
properly visualizing the accre tion flow and jet, we follow 
MacFadven & Woosley (1999) and show both the negative 
and positive x-region by duplicating the axisymmetric re- 
sult across the vertical axis. Color represents \og{po/ po.disk) 
with dark red highest and dark blue lowest. The final state 
has a density maximum of po w 2po,dM and a minimum of 
po ~ 10~ 13 po,£feA: at large radii. Grid zones are not smoothed 
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FIG. 3. — Jet has pummelled its way through presupernova core and through 
1/1 Oth of entire star. Time is t = 1.4 X 10 4 f g . Panel (A) shows final distribu- 
tion of logpo on the Cartesian plane. Black hole is located at center. Red is 
highest density and black is lowest. Panel (B) shows magnetic field overlayed 
on top of log of density. Outer scale is r = W^GM/c 2 . 



to show grid structure. Outer radial zones are large, but outer 
8 zones are below the resolution of the figure. 

Clearly the jet has pummelled its way through the surround- 
ing medium, which corresponds to the stellar envelope in the 
collapsar model. By the end of the simulation, the field has 
been self-consistently launched in to the funnel region and 
has a regular geometry there. In the disk and at the surface 
of the disk the field is curved on the scale of the disk scale 
height. Within r < \0 2 r g the funnel field is ordered and sta- 
ble due to the poloidal field dominance. However, beyond 
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FIG. 4. — Strongest magnetic field near black hole in X-configuration due 
to Blandford-Znajek effect and collimation of disk+coronal outflow. As in 
figure^] but outer scale is r= W 2 GM/c 2 . Time is t = 1.4 X 10 4 ( s . Black hole 
is black circle at center. Color scale is same as in figurel3l 



r ~ 10— 10 r g the poloidal field is relatively weak compared 
to the toroidal field and the field lines bend and oscillate er- 
ratically due to pinch instabilities. The radial scale of the 
oscillations is 10 2 r g (but up to 10 3 r g and as small as 10r ? ), 
where r ~ 10 r g is the radius where poloidal and toroidal field 
strengths are equal. By the end of the simulation, the jet 
has only fully evolved to a state independent of the initial 
conditions at r w 5 x 10 3 r g , beyond which the jet features 
are a result of the tail-end of the initial launch of the field. 
The head of the jet has passed beyond the outer boundary of 
r = I0 4 r g . Notice that the magnetic field near the black hole 
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R c 2 /(GM) 

FIG. 5. — Contours for t fa 1500f g and an outer scale of r Rj 10 2 r s , where 
the disk-corona boundary is a cyan contour where j3 fa «/fo 2 = 3, the corona- 
wind boundary is a magenta contour where (3 = 1, and the jet-wind boundary 
is a red contour where b 2 /(poc 2 ) = 1. Black contour denotes boundary be- 
yond which material is unbound and outbound (wind + jet). Beyond r fa I0r g 
the jet undergoes poloidal oscillations due to a toroidal pinch instability. 

is in an X-configuration. This is due to the BZ-effect having 
power Pj et oc sin 2 ^, which vanishes at the polar axis. The X- 
configuration is also related to the fact that the disk+corona is 
collimating the Poynting-dominated jet. The field is mostly 
monopolar near the black hole, and such field geometries 
decollimate for rapidly rotating black holes in force-free elec- 
trodynamics ( Krasnonolskv et al. 2005). 

Figures |5] and [6] show the energy structure of the disk, 
corona, a nd jet. T his figu re is comparab le to the left panel of 
figure 2 in McKinnev & Gammie (2004). Figure|5]shows con- 
tours for t « 1500?!, and an outer scale of r = I0 2 r g , while fig- 
ure|S]shows contours for (« 1 .4 x I0 4 t g and an outer scale of 
r = 10 3 r g . The disk-coronal boundary is represented as a cyan 
contour where j3 = p g /pb = 2(7- \)u/b 2 = 3, the coronal- 
wind boundary as a magenta contour where (3 = 1, and the 
jet- wind boundary as a red contour where b 2 /(p {) c 2 ) = 1. The 
black contour denotes the boundary beyond where material is 
unbound and outbound (wind or jet). Beyond r w 10r ? the jet 
undergoes poloidal oscillations due to toroidal pinch instabil- 
ities, which subside by r m 10 2 - 10 3 r ? where the jet is ther- 
mally supported. At large scales, the cyan and magenta con- 
tours closer to the equatorial plane are not expected to cleanly 
distinguish any particular structure. 

Figure Q shows the disk, corona, and jet magnetic field 
structure during the turbulent phase of accretion at t w 1500^. 
Compared to figure |4] this shows the turbulence in the disk, 
but is otherwise similar. The jet, disk, and coronal structures 
remains mostly unchanged at late times despite the decay of 
disk turbulence. That is, the current sheets in the disk do 
not decay and continue to support the field around the black 
hole that leads to the Blandford-Znajek effect. The corona 
thickness and radial extent do not require the disk turbulence, 
which only adds to the time-dependence of the coronal out- 
flow. 

Figure|8]shows a color plot of T^, where red is highest and 




R c 2 /(GM) 

FIG. 6. — Same as figurel5lbut for t fa 1 .4 X 10 4 ? g and an outer scale of r fa 
10 3 r s . Disk wind leads to broad coronal outflow. Jet remains collimated to a 
fixed opening angle. By r > 100r g , pinch instabilities subside once magnetic 
energy converted into thermal energy and supports jet. Residual slow dense 
and fast diffuse patches from the instability are present in the jet, such as the 
slower and cooler dense blob shown at top left corner of figure. 




r j0 SO 

R cVGM 

FIG. 7. — Field geometry near black hole for ; f» 1500< g during phase of 
strong disk turbulence. 

reaches up to r M - 10 3 - 10 4 and yellow has ^ - 10 2 - 10 3 . 
Panel (A) shows outer scale of r = I0 4 r g , while panel (B) 
shows outer scale of r = 10 3 r ? with same color scale. The 
inner-radial region is not shown since is divergent near 
injection region where ideal-MHD breaks down. Different re- 
alizations (random seed of perturbations in disk) lead to up 
to about Too ~ 10 4 as shown for the lower pole in the color 
figures. This particular model was chosen for presentation for 
its diversity between the two polar axes. The upper polar axis 
is fairly well-structured, while the lower polar axis has under- 
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FIG. 8. — Jet becomes conical at large radii with a core half-opening angle 
9j rj 5°. Plot of logToo with red highest (T^ ~ 10 4 ) and blue lowest. 
Yellow is ~ 10 2 - 10 3 . Panel (A) has outer scale of r = W 4 GM/c 2 . 
Panel (B) outer scale r = 10 3 GM/c 2 . Jet is independent of initial conditions 
by rzi5 X 10 3 r g . 

gone an atypically strong magnetic pinch instability. Various 
realizations show that the upper polar axis behavior is typi- 
cal, so this is studied in detail below. The strong hollow-cone 
structure of the lower jet is due to the strongest field being lo- 
cated at the interface between the jet and envelope, and this is 
related to the fact that the BZ-flux is oc sin 2 0, which vanishes 
identically along the polar axis. It is only the disk+corona that 
has truncated the energy extracted, otherwise the peak power 
would be at the equator. 
Figure [9] shows the velocity structure of the Poynting- 




1 10 10 2 10 3 10 4 

r c 2 /GM 



FIG. 9. — A radial cross section along a mid-level field of the jet showing 
the velocity structure. Top panel shows formation of magnetic fireball where 
matter (MA) and electromagnetic (EM) energy flux per unit rest-mass flux 
has ~ r V oo ; . 

dominated jet along a mid-level field line. The top panel 
shows the late-time time-averaged values of rJ^ M) (solid line), 
rj^ A) (dotted line), and Y = u'[BL- coords] (dashed line) as a 
function of radius along a mid-level field line. For 2 < r < 10, 
the value of r^ M) is highly oscillatory. This shows the lo- 
cation of the stagnation surface, where the poloidal veloc- 
ity u p = 0, is temporally variable. The stagnation surface 
varies between 2 < r < 10, within the range studied in prior 
sections. This is where ideal-MHD approximation is break- 
ing and thus at any moment the value of r^ M) nearly di- 
verges. Notice that while at r < 10 3 r ? the Poynting flux dom- 
inates, the Poynting flux energy is slowly converted into en- 
thalpy flux due to nonlinear time-dependent shock heating. 
Sho cks are expected beyond the m agneto-fast surface (see, 
e.g., Bogovalov & Tsinganos 2005) and here are due to pinch 
instabilities iEichledll993t lBegelmani n'998). For r > I0 3 r g , 
the enthalpy flux and Poynting flux are in equipartition. 

Thus a "magnetic fireball" has developed and the terminal 
Lorentz factor is of order ~ 400 for this choice of flow 
line. The Lorentz factor by r ~ 10 4 is still only 5 < Y < 10, 
with smaller patches with Y ~ 25. This implies there will be 
an extended acceleration region where the magnetic fireball 
loses energy during adiabatic expansion. These results should 
also he lp motivate the boundary conditions used in jet s imula- 
tions JAlov et alJ2000tlZhang. Wooslev. & M acFadve nl2003t 
Zhang W. et al.12 004). which sometimes start out with r ~ 50 
and u/po ~ 150 near the black hole. Such a high enthalpy 
only occurs by 10 3 r ? in our simulations. Their simulations are 
somewhat indicative of the Lorentz factor growth beyond our 
simulated time and radial range. This suggests that Y ~ 100 
should occur by r ~ 10 6 r f and that Y ~ 10 3 should occur by 
r ~ 10 7 r g . This is sufficient to avoid the compactness prob- 
lem. A simulation of the acceleration region is left for future 
work. 

The next lower panel of figure [5] shows the electromag- 
netic contribution to the specific angular momentum at infin- 
ity (L = QBj,) (solid line), the matter contribution to the spe- 
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FIG. 10. — A 9 cross-section of the jet at r = 5 X 10 3 GM/c 2 showing 
the velocity, density, energy, and magnetic structure in Gauss. Velocity and 
energy structure show Gaussian fits as dashed lines. Vertical dashed lines 
marks outer and core jet. 

cific angular momentum at infinity (hu^) (dotted line), and the 
specific angular momentum (dashed line). 

The bottom panel of figure [9] shows the orthonormal ra- 
dial (u r = y/g^u') (solid line) and cf> (u^ = ^Y<M> U ^) (dotted 
line) 4-velocities in Boyer-Lindquist coordinates. The di- 
rected motion becomes relativistic, while the ^-component 
of the 4-velocity becomes sub-relativistic. However, the ro- 
tation is still large enough that it may stabilize the jet against 
m = 1 kink instabilities, as shown in the nonrelativistic case 
{Nakamura & Meierl2004 . 

Figure [H)] shows the 9 cross-section of the jet at r = 5 x 
10 3 r ? for the upper polar axis according to figures[3] 0] and[8] 
This is a location by which the jet has stabilized in time, where 
farther regions are still dependent on the initial conditions. 
The hot fast "core" of the jet includes 9 < 5° at this radius 
and is marked by the left dashed vertical line. Also, it is useful 
to notice that 9 ss 0.14 corresponds to the radius obtained for 
the "mid-level" field line shown for radial cross-sections in 
figuresl9landll II 

The region within 9 < 27° is an expanded cold slow portion 
of the jet. It is possible that the gamma-ray burst photons are 
due to Compton drag from soft photons emitted by this jet 
"sheath" dumping into the faster spine (iBegelman & Sikora 
119871 iGhisellini et all2000t lLazzati et al.l2004l) . Based upon 
the data fits described below for r > I20r g and a radiation- 
dominated plasma, the outer sheath's (9 « 0.2) seed photon 
temperature as a function of radius is 

^~50kev(^) 1/3 (9) 

These seed photons can be upscattered by the jet and produce 
a GRB and high energy components. A self-consistent study 
of a Compton dragged jet that determines T and the emission 
energy is left for future work. 

The right dashed vertical line marks the boundary between 
the last field line that connects to the black hole. Beyond this 
region is the surrounding infalling medium. The upper-left 



panel shows (solid line), T (dotted line) and Gaussian fits 
to these two as dashed and long-dashed overlapping lines. The 
lower left panel shows the angular energy structure of the jet 
where e(9) = -r 2 T t r « ^poiiT^ with an overlapping Gaus- 
sian fit as a dotted line. The right panels show the density and 
magnetic structure of the jet. The upper right panel shows the 
density (solid line) and internal energy density (dotted line). 
The lower right panel shows the orthonormal Boyer-Lindquist 
coordinate toroidal (solid line) and radial (dotted line) field 
components. 

This jet structure is weakly due to an interaction with the 
surrounding medium, and primarily due to internal evolution 
of the jet. Clearly the slower jet envelope is non-negligible, 
giving c redence to the un iversal stru ctured jet (US J) mode l, 
but see lLamb et alJ ( 120041) and see IZhang B. et all ( 120041) : 
iLlovd- Ronning etall d2004l) . Gaussian jets have been used 
to expla in a universal connec tion between GRBs and x-ray 
flashes (Zha ng B. et alJ 2004). In most of our simulated 

models the jet structure is Gaussian with e{9) = eo<?~ e ^ w °, 
where erj « 0.18 and 9q ps 8°, within the range that they sug- 
gest fits observations. The total luminosity per pole is Lj « 
0.023Moc 2 , where 10% of that is in the "core" peak Lorentz 
factor region of the jet within a half-opening angle of 5°. 
Also, Too is approximately Gaussian with r^o «3x 10 3 and 
6*o ~ 4.3°. Also, r is approximately Gaussian with To ~ 5 and 
9i)SJ 11°. This can be folded into various models, such as the 
probability of observing polari sed emission in Compton dra g 
emission models (Ghisell mi et alJ2000t lLazzati et all2 004). 

For the collapsar model the above gives Li « 4 x 
10 51 ergs _1 for this j = 0.9375 black hole model. This is 
approximately equal to the luminosity emitted by the black 
hole plus the energy in pairs produced by annihilating neutri- 
nos. Not all of this energy can be observed i n 7-rays to ob- 
tain E j w2x 10 51 erg for a 30 second event (Zhang B. et alJ 
2004), unless the black hole has spin j sa 0.6 for most of the 
burst, which would s uggest little spin evolution that requires 
a much thicker disk ( Gammie, Shapiro, & McK innevll2004l) 
and so less neutrin o emission and annihilation efficiency. In 
McK innevI J2005b). we found that is weakly dependent 
on the black hole spin for a > 0.6, so the Lorentz factor is 
likely unaffected by such changes. 

It is more probable that the emission in 7-rays is ineffi- 
cient or only a fraction of the total energy is observed (i.e. 
this requires < 2% for each if sole effect) ; or the true black 
hole mass accretion rate is lower than predicted by current 
collapsar models (i.e. Mo ~ 0.002M Q .s~ 1 rather than Mo ~ 
O.IM0S" 1 if sole effect). For the standard collapsar model, if 
one only observed the core of the jet, which has only 10% of 
the luminosity, then an efficiency of converting to 7-rays of 
ss 20% is required to fit the model of Zhang B. etalJ (|2004). 

Notice that the disk thickness (and so mass accretion rate) 
may strongly determine the collimation angle. In higher 
mass accretion rate systems, the disk is thicker than described 
here, and the final opening angle may be smaller. Also, 
if the mass accretion rate is much lower, the disk is also 
thicker llKohri et alJ 120051) . The collapsar type model gives 
the thinnest disk near the black hole and small changes in 
the mass accretion rate lead to a small change in the disk 
thickness (Kohri et afl 120051) . Also, for the relevant black 
hole spins, th e spin only weakly determin es the final colli- 
mation angle (McKi nnev & Gammiel2004l) . The universal jet 
model requires quite large opening angles that would not be 
supported by invoking stellar (and so black hole) spin depen- 
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FIG. 1 1 . — A radial cross-section of both poles of (solid and dotted lines) 
jet along a mid-level field line in the jet showing the jet collimation, density, 

and magnetic structure (B^- r in Gauss). Overlapping dashed lines are fits. 

dence. It would also not be supported by invoking mass ac- 
cretion rate dependence, since a system with a hig her or lower 
accretion rate actually has a thicker disk ( Kohri et alJl2005l) 
and so stronger collimation. 

There is a weak dependence on the density of the 
stellar envelope compared to that s een in re lativis- 
tic hydrodynamic simulations of jets llAlov et al. J2| 
IZhang. Wooslev. & MacFadvedl2003t IZhang W. et alll200- 



since here the confinement is mostly magnetic. During this 
energetic phase of the jet, the medium plays little role in shap- 
ing the core of the je t, as has been suggested for other reasons 
jLazzati et alJ20 05). 

The peak Lorentz factors are obtained in a narrow region 
within a core with half-opening angle of 8 ss 5°. Notice 
as well that the very inner core of the jet is denser and 
slower, a f eature p redicted in the theoretical discussion of 
iMcKinnevI lH)05b). This feature is simply a result that the 
Blandford-Znajek flux is oc sin 2 #. Most of the Poynting en- 
ergy is initially at the interface between the jet and surround- 
ing medium, until the jet internally evolves to collimate into 
a narrower inner core. Notice that in McKinnev (2005g), we 
predicted ~ 1000 for ;" = 0.9, where the simulated results 
show a peak r m ps 3 x 10* and typical 10 2 <r m < 10 3 . This 
is in reasonable agreement with analytic estimates given the 
jet structure was not taken into account. 

Figure shows the collimation angle along a mid-level 
field line in the Poynting-dominated jet, along the same field 
line as in figure|9] and for both poles of the simulation. Gen- 
erally the two poles are similar but not necessarily identical. 
Many simulations that were performed show at least one polar 
jet develop some pinch instabilities. In the model shown, one 
side develops very strong instabilities by the end of the simu- 
lation. The development of the pinch instability mostly affects 
the far-radial Lorentz factor (T and — and so densities). 
The value of T is up to 5 times larger in pinched regions than 
non-pinched regions, while Too can be up to a factor 10 times 
larger in pinched regions compared to non-pinched regions. 

In figure for 7r ? < r < 100r g there is a region of col- 



limation slightly faster than logarithmic. For the field line 
closest to the coronal wind (not shown) starting at 9j w 1 .0, 
collimation is logarithmic with 



1 



2.8r c 



-1/3 



(10) 



Closer to the polar axis the collimation is faster. For the field 
line starting at 9 j ss 0.3, approximately 



-2/5 



(ID 



up to r ~ 10 2 r^. The inner-radial collimation is due to con- 
finement by the disk+corona, while in the outer-radial range 
the coronal outflow collimates the Poynting flow. The corona 
is likely required for the Poynting-dominated jet collima- 
tion, since without a disk, a monopole-like field near the 
black hole actually decol limates at higher black hole spin 
( Krasnopolskv et al. 2005). Far from the coronal outflow, the 
internal jet is collimated by internal hoop stresses. Note that 
the classical hoop-stress paradigm that jets can self-collimate 
is not fully tested here, but these results suggest that collima- 
tion is in large part due to, or at least requires, the coronal 
outflow. 

Notice from figure [9] that there is little acceleration be- 
yond r ~ 10 2 r f while the magnetic fireball is forming. It is 
well known that magnetic accele ration requires field lines to 
collimate (Begelman & Lil[l994l) . Once the flow is pseudo- 
conical the flow is unable to transfer magnetic energy into 
kinetic energy and proceeds to convert it into enthalpy flux. 
After r > 10 r g the flow oscillates around a pseudo-conical 
asymptote with a mean half-opening angle of 6j ~ 5°. The 
opening angle h as been found to be weakly d ependent on the 
black hole spin (McKinnev & Gammie 2004) and, as found 
in this paper, the details of the injection model. This resulting 
pseudo-conical asymptote results once the magnetic energy 
no longer dominates and the coronal wind no longer helps 
collimate the flow. This opening angle is in line with expecta- 
tions built around afterglow achromatic break measurements 
and estimates of the opening angle based upon energy of the 
burst. During the evolution of the burst, the jet may continue 
to collimate due to a more extended coronal outflow. This 
would be observed as an overall hard-to-soft evolution of 7- 
rays since one is gradually exposed to the slower edges of the 
jet since one likely does not observe exactly emission along 
the core center. The additional loading of neutrons from the 
coronal outflow may also lead to a hard-to-soft evolution of 
the jet as the overall Lorentz factor decreases over the dura- 
tion of the burst. 

Figure ^2 shows the toroidal field and pitch angle, which 
shows that eventually the toroidal field completely domi- 
nates the poloidal field, and the toroidal field remains or- 
dered. Thus, large po larizations up to 60% are possible 
(Granot & Konigl 2003). The inner radial toroidal field is 
well fit by 



= [Gauss] =0.0023 



( r 



The outer radial field is well fit by 



= [Gauss] =0.0023 



-0.7 



-1.5 



(12) 



(13) 



\J P0,diskC 2 ' 

The transition radius is r ss 390r g . Notice that this is along 
one mid-level field line, and the coefficients and power laws 
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are slightly different for each field line. The pitch angle shows 
that the field becomes toroidally dominated and at large ra- 
dius the magnetic loops have a pitch angle of < 1° by 10 3 r s . 
For regions not on the axis there is no p ossibility of recon- 
nection un less the flow is kink unstable (Drenkhahn 2002; 
iDrenkhahn & Sp ruit 2002). Given the regularity of the flow 
field and the lack of randomized field, reconnection is not nec- 
essary or likely. Here, the pinch instability drives a nonlinear 
coupling such as shocks that converts the Poynting flux to en- 
thalpy flux. This mechanism should be investigated in future 
work. 

The toroidal field dominance leads to m = pinch mag- 
netic instabilities. The m = 1 kink mode may also operate, but 
studying the 3D jet structure is left for future work. As such 
oscillations develop, this leads to arbitrarily sized patches 
moving at arbitrary relative velocities as in the internal shock 
model. The typical size of a patch is ~ 100 r g to ~ 10 3 r g in the 
lab frame along the length of the jet. For long-duration GRBs 
lasting about 30 seconds, the number of pulses should be no 
larger than the ratio of the scale of the instability to the length 
of the event (10 6 r ? ), or about 100- 1000 pulses. Different 
realizations of the same simulation show that the likelihood 
of an instability growing is random, which could lead to the 
large variations in the number of observed pulses and e xplain 
the diversity of observed pulses JNakar & Piran 2002). That 
is, some jets may have very few patches and so pulses. This 
type of instability is likely required to produce the diversity 
of GRBs. It is uncertain whether toroidal field instability in- 
duced variability dominates pair creation loading variability 
due to disk structure variability, which can be addressed by 
future work that self-consistently treats the neutrino emission 
from the disk. 

Figure^2 a l so shows the rest-mass density per disk density, 
internal energy density per rest-mass density, and the ratio of 
(twice) the comoving magnetic energy density per unit rest- 
m ass density. Notice that the density is close to that estimated 
in McKinnev (2005b), so those calculations are likely reason- 
able approximations to the simulation results. The inner radial 
rest-mass density is well fit by 
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The outer radial rest-mass density is well fit by 
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PO.disk 
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The transition radius is r 
internal energy density is moderately fit by 

-l. 



120r g . Likewise, the inner radial 



= 4.5 x 10" 



(inner) (16) 



P0,diskC* V 120r g 

The outer radial internal energy density is moderately fit by 

-1.3 



U 



= 4.5 x 10" 



(outer) (17) 



po,diskC" \ 120r g 

The transition radius is r « 120^ and is due to the presence 
of the thick disk at small radii. Notice that this is along one 
mid-level field line, and the coefficients and power laws are 
slightly different for each field line. One can compare this to 
the "envelope" density model used with 

-1.5 



PO 
PO.disk 



= 8x10" 



120r„ 



(envelope) (18) 




FIG. 12. — Jet becomes superfast and toroidal field dominates at large 
radii. Stagnation surface time-dependent but stable. Poynting-dominated jet 
characteristic (and other interesting) surfaces. Red lines are field lines. From 
r = outwards: Blue#l: horizon + ingoing-fast ; Cyan#l: ingoing-Alfven 

Black#l: B^ =B f ; Green#l: ergosphere ; Purple#l: ingoing-slow ; Black#2 
stagnation surface where poloidal velocity u 1 ' = ; Purple#2: outgoing-slow 

Cyan#2: outgoing- Alfven; Black#3: B$ =B' again ; Black#4: light cylinder 
; Blue#2: outgoing-fast. 



The enve lope has little impact on the jet structure. 

Figure fT2l shows the characteristic structure of the jet. The 
figure shows a log-log plot of one hemisphere of the time- 
averaged flow at late time. In such a log-log plot, 45° lines 
correspond to lines of constant 9. Lines of constant spherical 
polar r are horizontal near the z-axis and vertical near the R- 
axis. The field lines are shown as red lines. The disk and 
coronal regions have been truncated with a power-law cutoff 
for r < 100r g and a conical cutoff for larger radii. The inner- 
most blue line is the event horizon. Clearly the field lines are 
nearly logarithmic until r ~ 100r g . After this point the field 
lines stretch out and oscillate around a conical asymptote. 

Blue lines in figure [21 show the ingoing and outgoing fast 
surfaces. Clearly the transition to a supercritical (superfast) 
flow has occurred. Within r < 10r g the plotting cutoff creates 
the appearance that the fast surface and other lines terminate 
along the cutoff. 

The inner magenta line is the ingoing slow surface, while 
the outer magenta line is the outgoing slow surface. The in- 
ner cyan line is the ingoing Alfven surface. Energy can be 
extracted from the black h ole if and only if the Al fven surface 
lies inside the ergosphere (Takahashi et al. 1990). The ergo- 
sphere is shown by the green line, which is indeed outside the 
ingoing Alfven surface. The outer cyan line is the outgoing 
Alfven surface. The inner black line (next to the ergosphere) 

is the surface where = B r in Boyer-Lindquist coordinates. 
At larger radii there are two overlapping black lines. One 
corresponds to, again, where the toroidal field is equal to the 
poloidal field. The other corresponds to the light "cylinder" 
surface. 

Other small artifacts in the plot are a result of unsteady na- 
ture of the flow. However, the despite the unsteady nature of 
the flow the characteristic structure is quite simple and rela- 
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tively smooth in appearance. This indicates that the flow is 
mostly stationary. The stagnation surface is fairly unsteady, 
but stable. 

All the results discussed in this paper are robust to increases 
or decreases in numerical resolution. Convergence testing has 
been performed by choosing different models of the grid, 
including different Qj e/ parameters that controls whether the 
disk or jet is more resolved at small radii. Another resolu- 
tion of 256 x 256 instead of 512 x 256 was also chosen. No 
qualitative differences were found. The interpolation scheme 
was also changed from parabolic to linear and no qualitative 
differences for the far-field jet region were found. 

We find that the pair injection model weakly determines the 
bulk structure of the flow at large radii. Other injection mod- 
els were simulated by varying the total energy injected e', n/ , 
fraction of rest-mass created /„, fraction of internal energy 
created /},, fraction given to momentum energy f m , and the 
0-velocity of the injected particles. Only e',„ ; and f p strongly 
determine the bulk jet flow. 

Most of these results are insensitive to the two different 
models of the surrounding medium. One model is a surround- 
ing infall of material and the other is an "evacuated" exterior 
region. The jet structure is negligibly broader or narrower in 
the evacuated case due to magnetic confinement and collima- 
tion by the coronal outflow. 

3.6. Simulation as Applied to AGN and X-ray Binaries 

Notice that, quantitatively, all of these results can be simply 
rescaled by density in the jet in order to apply to other GRB 
models, AGN, and even approximately x-ray binary systems. 
This is because by the outer regions simulation, the jet has 
only reached about V ~ 5 - 10 and Too oc 1/ poj e t ■ 

In particular, numerical models for M87 give the same re- 
sults with 2 < r < 10 patches by r ~ 10 3 r g . The key dif- 
ference is that while the collapsar jet is very optically thick, 
the M87 jet is not so most of the shock-heated energy is lost 
to synchrotron radiation. The same structured jet is formed, 
and is useful to explain Te V BL Lac objects and radio galax- 
ies dGhisellin i et all2 005l) and may help explain observations 
of blazars ( Blandford & Levinsonlll995tlGhisellini & M adau 
U996l fChiaberg e et all2000l) . The full opening angle of the jet 
core of ~ 1 0° also agrees with the observations of the far- field 
jet in M87 (Junoret al. 1999). As discussed previously, the 
shocks are synchrotron cooled and the heat generated does not 
contribute to larger scale acceleration. While self-consistent 
synchrotron cooling should be included, these results suggest 
the shock zone (or "blazar zone" for blazars) should be quite 
extended between r ~ 10 2 r ? and up to about r ~ 10 r g . 

Similarly, simulations were performed using an injection 
model for GRS 1915+105 giving patches with 1.5 < T < 5 
by r ~ 10 2 r g . Due to the efficient pair-loading there is lit- 
tle extra Poynting flux to convert to enthalpy flux or kinetic 
energy fl ux by the time shoc ks occur at r ~ 10 2 r ? . As dis- 
cussed in McKinnev (2005b), the Poynting-lepton jet in the 
GRS 1915+105 model is possibly destroyed by Compton drag 
by disk photons and synchrotron self-Compton and limited to 
at most T < 2, but the jet could be completely destroyed by 
Compton drag. 

3.7. Summary of Fits 

A summary of the fits along a fiducial field line is given. 
Near the black hole the half-opening angle of the full 
Poynting-dominatedjet is 6j ~ 1.0, while by r ~ \20r g , Oj ~ 



0.1. This can be roughly fit by 



-0.4 



(inner) 



(19) 



for r < I20r g and Oj ~ 0.14 beyond. The core of the jet fol- 
lows a slightly stronger collimation with 

6j oc r" 2/5 (20) 

up to r < 120r g and Oj — 0.09 beyond. Also, roughly for M87 
and the collapsar model, the core of the jet has 

/ r \ 0.44 

— J (inner) (21) 

for 5 < r < 10 3 r g and constant beyond for the M87 model 
if including synchrotron radiation, while the collapsar model 
should continue accelerating and the power law will trun- 
cate when most of the internal and Poynting energy is lost 
to kinetic energy and the jet becomes optically thin at about 
r ~ 10 9 r g or internal shocks take the energy away. If the 
acceleratio n is purely thermal with out any magnetic effect, 
then r«r (M eszaros & Reeslll997l) . However, it is not clear 
how the equipartition magnetic field affects the acceleration. 
Roughly for GRS 1915+105 the core of the jet has 

( r \° H 

r bulk ~ (inner) (22) 

for 5 < r ~ 10Vg and constant beyond, with no account for 
Compton drag or pair annihilation. Also, for any jet system 
the base of the jet has po oc r" 0,9 (inner) for r < \20r g and 
oc r" 2,2 (outer) beyond. For the collapsar and M87 models 
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(outer), 
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P0,disk \ xi - u '«, 

while for GRS 1915+105 the inner-radial coefficient is 10~ 5 
and outer is 6 x 10~ 3 . For the collapsar model, the inner radial 
internal energy density is moderately fit by 
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PO,diskC* 



= 4.5 x 10" 
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(inner). (25) 



The outer radial internal energy density is moderately fit by 
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= 4.5 x 10" 9 



(outer). (26) 



Pa,diskC' \ 120r g 

The transition radius is r « \20r g . For M87 the internal en- 
ergy is near the rest-mass density times c 2 until r ~ \20r g 
when the dependence is as for the collapsar case. For 
GRS 19 15+ 105 the internal energy is near the rest-mass den- 
sity times c 2 until r ~ \20r g and then rises to about 2.5 times 
the rest-mass density times c 2 . The inner radial toroidal lab 
field is well fit by 



V P0,diskC 



: [Gauss] =0.0023 



390r e 
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(inner) (27) 



for 5 < r < 390r g . The outer radial toroidal lab field is well fit 



by 
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for r > 390r g . 

For the typical jet with no atypical pinch instabilities, the 
energy and velocity structure of the jet follow 

e(6) = e e- el l w K (29) 

where eo sa 0.18 and Oq sa 8°. The total luminosity per pole 
is Lj fts 0.023M c 2 , where 10% of that is in the "core" peak 
Lorentz factor region of the jet within a half-opening angle of 
5°. Also, Too is approximately Gaussian 

T oc (9) = T oofi e- e2 t w K (30) 

where o ~3x 10 3 and 9q w 4.3°. Also, T is approximately 
Gaussian 

T(0) = Y Q e- 'l w \ (31) 

where To ~ 5 and 6q ss 1 1°. The outer sheath's (9 ss 0.2) seed 
photon temperature as a function of radius is 

r 7 ,,^50kev(^-) 1/3 . (32) 

4. DISCUSSION 

A collapsar-type GRMHD simulation with a neutrino an- 
nihilation model and Fick diffusion model has been studied. 
A self-consistent Poynting-dominated jet is produced and the 
Lorentz factor at large distances is ~ 10 3 in the core of 
the jet, but an equally important structured component exists 
with Too ~ 10 2 . The Lorentz factor of the jet is determined by 
the electron-proton loading by Fick diffusion of neutrons. No- 
tice that estimates of the baryonic mass-loading from the op- 
tical flash of GRB990123 sug gest a baryon loading that gives 
r ~ 10 3 (Soderberg & Ramirez-Ruiz 2003), which is in basic 
agreement with the findings here. The half opening angle of 
the core of the jet is 9j ~ 5°, while there exists a significant 
structured component with a half opening angle of 9j < 25°. 
It is likely that this jet component survives traversing the re- 
maining portion of the stellar envelope (r ~ 10 5 r g , another 
decade in radius). 

The jet at large distances is unstable to, at least, an m = 
pinch instability due to the do minance of the t oroidal field 
to the poloidal field (see, e.g. Begelman 1998). Here it is 
found that the energy of the jet is patchy. Typical realizations 
of the same model have a random number of patches with 
100 < Too < 10 3 . This likely results in a pulsed sequence 
of magnetic fireballs that move at large rela tive Lorentz fac- 
tors, as required by the internal sho ck model (Meszaros 2002; 
iGhirlanda et alJl2003t lPiranll2005h . The acceleration region 
and internal shock region was not simulated since the dynam- 
ical range required is another ~ 6 orders of magnitude in ra- 
dius. 

Some models assume the Poynting energy to be con- 
verted into radiation fa r from the collaps i ng star by in ternal 
dissip ati on (see, e.g.. iTh ompsonl 119941: IMeszaros & Reesl 
19971 ISpruit. Daigne. & Drenkhahnl 1200 ll iDren khahn 
20021 Drenkhahn & Spruit! 120021 ISikoraetalJ 120031 
Lvutikov, Pariev, & Blandford 2003) such as magnetic 
reconnection. However, in our model, significant internal 
dissipation occurs already by 10 3 r g . 

A very efficient model that reduces the compactness prob- 
lem invokes Compton drag (bulk Comptonization) to gen- 
erate the emission, where the stellar envelope or presuper- 
nova region provi des soft photons (Ghis ellini et al.l l2000t 
lLazzati et alJl200 4). Our numerical model shows a jet struc- 
ture with a slow cold "sheath" necessary for supplying the 



seed photons, so this process might explain the GRB emis- 
sion pro cess. 

As in iMcKinnev & Gammid (|2004), a mildly relativistic 
Poynting-baryon jet is launched from the inner edge of the 
disk with a half-opening angle of roughly 16° to 45°, how- 
ever long-time study of this jet component depends on sus- 
taining long-time disk turbulence, which cannot be simulated 
in axisymmetry. For GRBs, this hot baryon-loaded jet might 
play a significant role in the subsequent supernova (by pro- 
ducing , e.g., 56 Ni) and be an interesting site for the r-process 
(Kohri et all 120051) . The mass ejection rate is found to be 
MQ tCorona i « 0.03Mo, which for a 30 second event gives 0. 1M Q 
of baryonic mass. This may be sufficient to power super- 
novae. The Lorentz factor is about Y ~ 1.5-3, such as o b- 
served in a component of SN1998bw (Kulkarni et al. 1998). 

For GRBs, the picture that emerges is that Poynting flux 
is emitted from the black hole at r ~ r g and is loaded with 
electron-positron pairs within r < 20r g . A similar mass- 
fraction of neutrons Fick-diffuses into the jet, and they dom- 
inate the rest-mass once the temperature decreases to T < 
6 x 10 9 K by r ~ 10r g . Magnetic acceleration occurs over 
r < 10 2 r„ up to r ~ 10. Notice that this is in contrast 
to IMeszaros & Reesl i 19971) who assumed the initial bulk 
Lorentz factor is Y[r = rg] ~ r ( ^ M) , which assumed instan- 
taneous magnetic acceleration. 

Once the toroidal field dominates the poloidal field, mag- 
netic instabilities develops around r ~ 10- 10 2 r g which turns 
the jet into an equipartition "magnetic fireball" by r ~ 10 3 r g . 
The spatial structure of the magnetic fireball energy flux is 
patchy on the instability scale of r ~ 10 2 r g . Simulations that 
demonstrate how these features would later interact require 
much more radial dynamic range. Between I0 4 r g < r < 10 8 r g 
the flow accelerates due to adiabatic expansion with a radial 
dependence up to Y oc r (Meszaros & Rees 1997), where the 
Poynting flux provides a reservoir of magnetic energy that is 
continuously shock-converted into thermal energy. 

After passing the stellar surface at r ~ 10 5 r g , the plasma 
becomes transparent to 7-rays at 10 6 r g < r < 10 7 r g , where 
possibly Compton drag of the sheath seed photons produce s 
the GRB emission (Ghise llini et alJ2000l:lLazzati et alJ2004l) . 
By r ~ 10 7 , patches of varying Y ~ 100- 1000 have reached 
their terminal velocity. Between 10 8 r g < r < 10 10 r g the fire- 
ball is optically thin to Compton scattering, radiation escapes, 
and pairs no longer annih ilate but can con tinue to be magnet- 
ically accelerated (Meszaros & Rees 1997). Within this same 
radial range, internal shocks proceed to convert the remain- 
ing kinetic energy to nonthermal 7-rays. Beyond 10 1( V g a re- 
verse shock may occur and external shocks with the interstel- 
lar medium (ISM) occur. 

Simulations of AGN were performed, which show similar 
r vs. radius as the collapsar case. That is, there is an early 
transfer of Poynting flux to kinetic energy flux leading up to 
about r ^ 5- 10 by about r ~ 10 3 r g for an M87 model. This is 
consistent with the l ack of observed Co mptonization features 
in blazars (see, e.g., Sikora et al. 2005), although this is also 
consistent with the fact that the optical depth to Comptoniza- 
tion is low in such systems. 

While prior work assuming r^/j ~ 10 and 6j ~ 15° found 
that the B Z power is insufficient to account for Blazar emis- 
sion (Maraschi & Tavecchio 2003), as shown here the struc- 
ture of the Poynting-dominated jet is nontrivial. The re- 
gion with Ybuik ~ 10 is narrower with 9j ~ 5° and jet emis- 
sion is likely dominated by shock accelerated electrons with 
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r e > 100 with an extended high energy tail. This lowers the 
necessary energy budget of the jet to be consistent with BZ 
power driving the jet. 

4.1. Theoretical Contemporaneous Comparisons 

Theoretical studies of ideal MHD jets focused on cold 
jets and the conversio n of Poynting energy directly to ki- 
netic energy (see, e.g ., |LLe t alJll9 92: Begelman & Lil ll994l 
Daigne & Drenkhahn 2002). Indeed, much of the work 
has focused on self-similar models, of which the only self- 
consistent solutio n foun d is suggested to be R-self-similar 
models ( Vlahakisl l2004l) . Unfortunately, such models result 
in only cylindrical asymptotic solutions, which is apparently 
not what is observed in AGN jets nor present in the simula- 
tions discussed in this paper. The previous section showed 
that some aspects of the jet are nearly self-similar and that 
there is some classic ideal-MHD acceleration occurring in 
this region. However, once the toroidal field dominates the 
poloidal field, the Poynting flux is converted into internal en- 
ergy until they reach equipartition. The conversion is due to 
a time-dependent nonlinear coupling, such as shocks, due to 
the magnetic pinch instability. Cold ideal-MHD jet models 
would have no way of addressing this. This region also in- 
volves relatively rapid variations in the flow, so stationary jet 
models would have difficulty modelling this region. 

Notice that GRMHD numerical models show that the mag- 
netic field corresponding to the Blandford- Znajek effect dom- 
inates all other magnetic field geometries ( Hiros e~et all 2004; 
lMcKinnevl2005al) . For example, there are no dynamically sta- 
ble or important field lin es that t ie the black hole to the outer 
accretion disk as in lUzdenskvl J200l . All black hole field 
lines tie to large distances or tie to the horizon-crossing accre- 
tion disk. Also, there are no field lines that connect the inner- 
radial accretion disk and large distances. The Blandford- 
Znajek associated magnetic field is completely dominant, as 
shown in, e.g., figure[3]and figure^ 

We have not directly included t hose effects of neutron dif - 
fusion that induce a jet structure tLevinson & Eichlerll2003l) . 
However, there is probably strong mixing within the jet and 
the large-scale jet is probably not affected by the details of 
the interface where neutron-diffusion occurs. Alternatively, 
the Gaussian structure we find may dominate the power law 
structure they find since the Gaussian structure is due to a 
strong internal electromagnetic evolution of the jet and matter 
plays little role in setting the jet structure. 

4.2. Numerical Contemporaneous Comparisons 

MacFadve n & Wooslevl ( fT999i) evolve a presupernova core 
with a black hole inserted to replace part of the core and 
evolve the nonrelativistic viscous hydrodynamic equations of 
motion for various values of viscosity coefficient (a), mass 
accretion rate, nuclear burning, stellar angular momentum, 
and some models have an injected energy at the poles at some 
fixed energy rate. For models in which they inject energy 
at some specified rate, they find the baryon-contamination 
problem is somewhat alleviated by forming a hot bubble. As 
shown in their figure 28, they find the jet has u/(poc 2 ) ~ 10 
and they suggest that this corresponds to ^10, insuffi- 
cient to avoid the compactness problem. Their a = 0.1 model 
shows a significant disk outflow, which over the duration of 
the GRB would yield M ~ M© in mass that could power a 
supernova. In their a = 0.01 model there are insignificant out- 
flows. 



Our MHD results are most comparable to their a = 0.01 
model. We find a disk outflow yielding M ~ 0.1M Q , which 
may still be sufficient to produce a supernova component. We 
have self-consistently evolved the jet formation and included 
an approximate model for the pair creation physics. In our 
case the baryon-contamination problem is self-consistently 
avoided by magnetic confinement of the jet against baryons, 
where only a small number of neutrons Fick-diffuse into the 
jet and self-consistently determine the Lorentz factor at large 
distances. We find that neutrino annihilation energy is proba- 
bly dominated by energy from the black hole. Notice that we 
find Too ~ 100- 1000, which is much larger than they find. 
The difference between their and our model is the presence 
of a magnetic field and a rotating black hole, which drive the 
B Z-effect and a stronger evacuation of the polar jet region. 

Proga (2005) have performed nonrelativistic MHD simula- 
tions of collapsars with a realistic equation of state. They sug- 
gest MHD accretion is able to launch, collimate, and sustain a 
strong Poynting outflow, although they measure a jet velocity 
of v ~ 0.2c. They find the jet is Poynting flux dominated such 
that Too < 10, insufficient to avoid the compactness problem. 
The rotation of the black hole is crucial to generate a suffi- 
ciently Poy nting-dominated jet. S imilar fully relativistic sim- 
ulations by Mizuno et al. (2004) have performed with a jet 
velocity of only v ~ 0.3c, which is due to their choice of the 
"floor" model in the polar regions and the short time of inte- 
gration. 

IZhang. Wooslev. & MacFadvenl ( 120031) : IZhang W etai] 
(2004) use a relativistic hydrodynamic model to simulate 
jet propagation through the stellar envelope and subsequent 
breakout through the stellar surface. They assume a highly 
relativistic jet is formed early near the black hole and they 
inject the matter with a large enthalpy per baryon of about 
u/poc 2 ~ 150. Th ey also tune the in jected energy to compare 
with observations (iFrailet al.l2001l) . In their view, variability 
is due to hydrodynamic instabilities between the cocoon and 
core of the jet. They find the stellar envelope can collimate 
the flow. 

We have self-consistently evolved the formation of the jet 
along with the disk without having to assume a jet struc- 
ture. We suggest that the injected energy per baryon is only 
u I poc 2 < 20 unless super-efficient neutrino mechanisms are 
invoked. We also find that the self-consistent energy released 
is larger than observed, suggesting an fairly inefficient genera- 
tion of 7-rays and a dominant lower T jet component that may 
explain x-ray flashes. We suggest magnetic instabilities due 
to pinch or kink modes dominate the variability rather than 
hydrodynamic instabilities, which are known to be quenched 
by magnetic confinement. We find that the jet structure is 
negligibly broader for models with no stellar envelope due to 
magnetic confinement and we find that the outer portion of the 
Poynting-dominated jet is collimated by the coronal outflow. 

Many hydrodynamic jet propagation simula ti ons have been 
studied (for a review see iScheck et alJl2002l) . iScheck et alJ 
(2002) evolve hydrodynamic relativistic hadronic and lep- 
tonic jets using a realistic equation of state. They find that the 
resulting morphology of the jets are similar, despite the differ- 
ent composition. This suggests jet gross morphology is not a 
useful tool to differentiate jet composition. In our model, the 
most relativistic jet component is lepton-dominated in AGN 
and x-ray binaries and baryon-dominated in GRB systems. 
Notice that magnetic confinement quenches most hydrody- 
n amic instabilities. 

iDe Villiers et alJ (|2005b) evolve a fully relativistic, black 
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hole mass-invariant model, showing a jet with hot blobs mov- 
ing with r ~ 50. A primary result in agreement with the 
results here is that a patchy or pulsed "magnetic fireball" is 
produced. This is gratifying since it suggests that the devel- 
opment of a "magnetic fireball" is not an artifact of the nu- 
merical implementation but is a likely result of shock heating. 
We also agree in finding that the core of the jet is hot and fast 
and is surrounded by a cold slow flow. One difference is that 
they say they seem to find the flow is cylindrically collimated 
by t ft 300r g , while we find nearly logarithmic collimation 
until r ~ I0 2 r g and an oscillatory conical asymptote beyond. 
Another difference is that they suggest temporal variability is 
due to injection events near the black hole, while we suggest 
it is due to pinch (or perhaps kink) instabilities at r > I0 2 r g . 

Another difference is that they find a larger value of T at 
smaller radius than we find. Indeed, one should wonder why 
their black hole mass-invariant model applies to only GRBs 
and not AGN or x-ray binaries. This is a result of their much 
lower "floor" density of po ~ 10~ 12 po,disib which is not con- 
sistent with self-consistent pair-loading or neutron diffusion 
loading. Also, because they use a constant density floor, the 
jet region at large radius must be additionally loaded with an 
arbitrary amount of rest-mass. The typical "floor" model adds 
rest-mass into the comoving frame, but this artificially loads 
the jet with extra mass moving at high V. We suggest that 
rather than the acceleration occurring within r < 10Qr g of their 
simulation, that acceleration occurs much farther away before 
the emitting region at r < I0 9 r g . In their T ~ 50 knots the 
gas has 10 3 < u/poc 2 < 10 6 . This implies that their actual 
terminal Lorentz factor is on the order of 10 5 < < 10 7 , 
which would imply that the external shocks occur before in- 
te rnal shocks could o ccur. 

iKomissarovl |2005) study the BZ process and MHD Pen- 
rose process. They found that the p rior "jet" results of 
llKoide. Shibata. Kudoh. & Meierl20 02). who evolve only for 
t ~ 100f ? , are only a transient phenomena associated with 
the initial conditions. They also suggest that there are seri- 
ous problems for the MHD-driven model of jets since they 
do not find jets in their models except for contrived field ge- 
ometries. However, this is potentially due to three effects. 
First, their outer radius is r ~ 50r g , while here we study out to 
r ~ 10 4 r g . The magnetic acceleration only leads to a logarith- 
mic increase in the velocity, so a large radial range is required 
to observe relativistic motion. S econd, magnetic accel eration 
requires collimating field lines (Begelman & Li 1994), and a 
disk or disk w ind is probably requir ed to collimate the mag- 
netic outflow jOkamotoll 1 9991 120001) . The "disk" that forms 
in their models is relatively thin. Third, to avoid numerical 
errors, they limit their solution to b 2 /(poc 2 ) < 100. Here we 
find that a self-consistent source of jet matter from pair cre- 
ation and neutron diffusion leads to b 2 /(poc 2 ) ~ 10 5 near the 
black hole. These three effects are probably why they find no 
relativistic jets. 

4.3. Limitations 

As has been pointed out by IKomissarovl J2005I) . noncon- 
servative GRMHD schemes often overestimate the amount 
of thermal energy produced. All GRMHD numerical models 
suffer from some numerical error. Shock-conversion of mag- 
netic energy to thermal energy is modelled by HARM in the 
perfect magnetic fluid approximation with total energy con- 
served exactly. However, our numerical model may overesti- 
mate the amount of magnetic dissipation in shocks, and so the 



shocks may more slowly convert magnetic energy to thermal 
energy. If this dissipation was delayed until the 7-ray photo- 
sphere at r ~ 10 7 - 10 8 r s , then those shocks could be directly 
responsible for the 7-ray emission from GRBs. Clarification 
of this issue is left for future work. However, we expect that 
toroidal field instabilities drive efficient magnetic dissipation 
as shown in the numerical results presented here. 

In order to evolve for a longer time than simulated in this 
paper, other physics must be included. For long-term evolu- 
tion of a GRB model, one must include disk neutrino cool- 
ing, photodisintegration of nuclei, and a realistic equation of 
state. If one wishes to track nuclear species evolution, a nu- 
clear burning reactions network is required. For the neutrino 
optically thick region of the disk, radiative transport should 
be included. The self-gravity of the star should be included 
to evolve the core-collapse. This includes a numerical rel- 
ativity study of the collapse of a rotating magnetized mas- 
sive star into a black hole (see review by Stergioulas 2003, 
§4.3). The jet should be followed through t he entire star and 
beyond penetration of the stellar surface (|A lov et al J l2000t 
Zhan g. Wooslev. & M acFadve nl2003tlZhang W. et al.l 2004). 

As applied to all black hole accretion systems, some 
other limitations of the numerical models presented include 
the assumption of axisymmetry, ideal MHD, and a non- 
radiative gas. The assumption of axisymmetry is likely 
not import ant for the inner jet reg ion since our earlier re- 
sults ( McKinnev & Gammie 2004) find quantitative agree - 
ment with 3D results JDe Villiers. Hawlev. & Krolikl 2003). 
The primary observed lim itation of axis ymmetry appears to 
be the decay of turbulence (Cowling 1934), which we attempt 
to avoid by requiring a resolution that gives quasi-steady tur- 
bulence for much of the simulation. Also, the jet at large dis- 
tances has already formed by the time turbulence decays, and 
by that time the jet at large radius is not in causal contact with 
the disk. 

When the toroidal field dominates the poloidal field, even- 
tually m = 1 kink instabilities and higher modes may appear 
(see, e.g., fNakamur a et alJl200lt iNakamura & Meier! 2004). 
Thus our models may underestimate the amount of oscillation 
in the flow and the conversion of Poynting flux to enthalpy 
flux. 

The limitation of axisymmetry may also limit the efficiency 
of the Blandford-Znajek process. The magnetic arrested disk 
(MAD) model suggests that any accretion flow likely accumu- 
lates a large amount of magnetic flux near the black hole. This 
means the efficiency of extracting energy wou ld be higher 
( Igumenshchev et al. 2003; Narava n et al.l2003l) . 

We have neglected low-energy jet and disk radiative pro- 
cesses in the numerical simulations. Future work should in- 
clude a model of synchrotron radiation for an AGN model in 
or der to simulate the " jet" emission to verify the claims made 
in McKinnev (2005b) that the broad "jet" emission observed 
by lJunor et al] i 19991) is actually the coronal outflow rather 
than the Poynting-dominatedjet. 

Also for AGN, a self-consistent simulation with syn- 
chrotron emission would likely show the continuous loss of 
Poynting flux until the synchrotron cooling timescale is longer 
than the jet propagation timescale, which still suggests the jet 
Poynting flux is finally in equipartition with the enthalpy flux. 
This would suggest that the shock-zone, and so emission re- 
gion, is more extended (r ~ 10 2 - 10 r g ) than the simulated 
shock-zone (r ~ 10 2 - 10 3 r g ). 

For AGN and x-ray binaries, the radiatively inefficient disk 
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approximation, which assumes electrons couple weakly to 
ions, may not hold. If the electrons and ions eventually couple 
near the black hole, then the disk might collapse into an un- 
stable magnetically dominated accretion disk (MDAF) ( Meier 
2005). Like the MAD model, this might drastically alter the 
results here, although it is uncertain whether jets are actually 
produced under the conditions specified by the MDAF model. 

The single-fluid, ideal MHD approximation breaks down 
under various conditions, such as during the quiescent out- 
put of x-ray black hole binaries, where a two-temperature 
plasma likely forms near the black hole as ions and electrons 
decouple (see, e.g., Naravan & Yi 1995). Resistivity plays a 
role in current sheets where reconnection events may gener- 
ate flares as on the sun, such as possibly observed in Sgr A 
JGenzel et alj |2003). Finally, radiative effects may introduce 
dynamically important instabilities in the accretion disk (e.g. 
lGammiell998tlBlaes & Socratesl2003l) . 

5. CONCLUSIONS 

Primarily two types of relativistic jets form in black hole 
(and perhaps neutron star) systems. The Poynting-dominated 
jet region is composed of field lines that connect the rotating 
black hole to large distances. Since the ideal MHD approxi- 
mation holds very well, the only matter that can cross the field 
lines are neutral particles, such as neutrinos, photons, and free 
neut rons. 

In lMcKinnevl ll2005b). we showed that the primary differ- 
ences between GRBs, AGN, and black hole x-ray binaries 
is the pair-loading of the Poynting-dominated jet, a similar 
mass-loading by free neutrons in GRB-type systems, the opti- 
cal depth of the jet, and the synchrotron cooling timescale of 
the jet; 

In lMcKimievl j2005bl) . we showed that for GRB-type sys- 
tems the neutron diffusion flux is sufficiently large to be dy- 
namically important, but small enough to allow T ~ 100 — 
1000. Beyond r ~ 10r g many of the electron-positron pairs 
annihilate, so the Poynting-dominated jet is dominated in 
mass by electron-proton pairs from collision-induced neutron 
decay. Most of the energy is provided by the BZ effect instead 
of ne utrino-anni h ilation. 

In lMcKinnevI (12005b). we showed that for AGN and x- 
ray binaries, the density of electron-positron pairs established 
near the black hole primarily determines the Lorentz factor at 
large distances. Radiatively inefficient AGN, such as M87, 
achieve 2 < Too ;$ 10 and are synchrotron cooling limited. 
The lower the 7-ray radiative efficiency of the disk, the more 
energy per particle is available in the shock-zone. Radia- 
tively efficient systems such as GRS1915+105 likely have no 
Poynting-lepton jet due to strong pair-loading and destruction 
by Comptonization by the plentiful soft photons for x-ray bi- 
naries with optically thick jets. However, all these systems 
have a mildly relativistic, baryon-loadedjet when in the hard- 
low state when the disk is geometrically thick, which can ex- 
plain jets in most x-ray binary systems. 

A GRMHD code, HARM, with pair creation physics was 
used to evolve many black hole acc retion disk models. The 
basic theoretical predictions made in McKinnev (2005tJ that 
determine the Lorentz factor of the jet were numerically con- 
firmed. However, Poynting flux is not necessarily directly 
converted into kinetic energy, but rather Poynting flux is first 
converted into enthalpy flux into a "magnetic fireball" due to 
shock heating. Thus, at large distances the acceleration is pri- 
marily thermal, but most of that thermal energy is provided 
by shock-conversion of magnetic energy. In GRB systems 



this magnetic fireball leads to thermal acceleration over an 
extended radial range. The jets in AGN and x-ray binaries 
release this energy as synchrotron and inverse Compton emis- 
sion and so the jet undergoes negligible thermal acceleration 
beyond r ~ 10 2 -10 3 r g . 

Based upon prior theoretical ( M cKinnevll2005bT) and this 
numerical work, basic conclusions for collapsars include: 

1 . Black hole energy, not neutrino energy, typically pow- 
ers GRBs. 

2. Poynting-dominated jets are mostly loaded by e~e + 
pairs close to the black hole, and by e~p pairs for 
r > I0r g . 

3. BZ-power and neutron diffusion primarily determines 
Lorentz factor. 

4. Variability is due to toroidal field instabilities. 

5. Poynting flux is converted into enthalpy flux and leads 
to the formation of a "magnetic fireball." 

6. Patchy jet develops 10 2 < Too < 10 3 , as required by 
internal shock model. 

7. Random number of patches (< 1000 for 30 second 
burst) and so random number of pulses. 

8. Energy structure of jet is Gaussian with 6q w 8°. 

9. Core of jet with 0j 5° can explain GRBs. 

10. Extended slower jet component with 6j w 25° can ex- 
plain x-ray flashes. 

1 1 . Coronal outflows with T ~ 1 .5 may power supernovae 
(by producing, e.g., 56 Ni) with M ~ O.IMq processed 
by corona. 

Based upon prior theoretical ( M cKinnevll2005bl) and this 
numerical work, basic conclusions for AGN or x-ray binaries 
include: 

1. Poynting-dominated jets e~e + pair-loaded unless advect 
complicated field. 

2. 7-ray radiative efficiency, and so pair-loading, deter- 
mines maximum possible Lorentz factor. 

3. Poynting-leptonjet is collimated with 9j « 5°. 

4. Extended slow jet component with 6j < 25°. 

5. For fixed accretion rate, variability is due to toroidal 
field instabilities. 

6. Poynting flux is shock-converted into enthalpy flux. 

7. In some AGN, shock heat in transonic transition lost 
to synchrotron emission and limits achievable Lorentz 
factor to 2 < T < 10 (e.g. in M87). 

8. Coronal outflows produce broad inner-radial jet fea- 
tures in AGN together with well-collimated jet compo- 
nent (e.g. in M87). 

9. In some x-ray binaries, Compton drag loads Poynting- 
lepton jets and limits Poynting-leptonjet to T < 2 or jet 
destroyed. 
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10. In some x-ray binaries, Poynting-lepton jet optically 
thick and emits self-absorbed synchrotron. 

1 1 . Coronal outflows have collimated edge with T < 1.5. 

12. Coronal outflows may explain all mildly relativistic and 
nonrelativistic jets in radiatively efficient systems (most 
x-ray binaries). 

For AGN and X-ray binaries, the coronal outflow collimation 
angle is strongly determined by the disk thickness. The above 
conclusions regarding the collimation angle assumed H/R ~ 
0.2 near the black hole andH/R ~ 0.6 far from the black hole, 
while H/R ~ 0.9 (ADAF-like) is perhaps more appropriate 



for some systems. The sensitivity of these results to H/R is 
left for future work. 
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APPENDIX 

A. GRMHD EQUATIONS OF MOTION WITH PAIR CREATION 

A single-component GRMHD approximation that accounts for baryon conservation is summarized. Leptons are assumed to 
be conserved in the equation of state (EOS) and by accounting for free-streaming neutrinos. We assume an ideal gas EOS 
of relativistic particles and neglect radiative transport since we assume all particles either stream freely or are trapped in the 
fluid. Alternatively, we assume the initial conditions well-model the steady-state radiative equilibrium. We model the streaming 
photon/neutrino-annihilation into pairs by injecting rest-mass and energy-momentum at an appropriate fraction of the baryon 
density. 

The black hole has a Kerr metric written in Kerr-Schild coordinat es, where the Kerr metric in K err-Schild coordinates and 
the Jacobian transformation to Boyer-Lindquist coordinates is given in M cKinnev & Gammid (120041) . We use Kerr-Schild rather 
than Boyer-Lindquist because in Kerr-Schild coordinates the inner-radial boundary can be placed inside the horizon and so out 
of causal contact with the flow. It is difficult to avoid interactions between the numerical inner boundary and the jet when using 
Boyer-Lindquist coordinates. This interaction leads to excessive variability in the jet since the ingoing superfast transition is not 
on the grid and then the details of the boundary condition significantly impact the jet. Numerical models of viscous flows have 
historically had similar issues (see, e.g. JMcKinnev & Gam mie 2002|). 

A single-component MHD approximation is assumed such that baryon number is conserved up to a source term due to pair 
creation due to either radiative annihilation or neutron diffusion, such that 

(p u^ = S p , (Al) 

where po = mi,rib, nib ~ ni„ the neutron mass, n/, = n n +n p , and m p is the proton mass. One may choose an arbitrary mass weight 
to define po in the single component fluid approximation. 
For a magnetized plasma the conservation of energy-momentum equations are 



f T fj,u T nv\ 
v-'ma^-'emJ 



(A2) 



where r* 1 " is the stress-energy tensor, which can be split into a matter (MA) and electromagnetic (EM) part. In the fluid 
approximation 

TZ = (po + u g )^u» + p g P^, (A3) 

with a relativistic ideal gas pressure p g = (7- \)u g , 7 = 4/3, and P^ v = + u^u v is the projection tensor. Either 7 = 5/3 or 
7 = 4/3 for either the bary on-dominated component or lepton-dominated component does not change the results described in 
section[3] 

A.l. Injection Physics 

A detailed pair creation model would self-consistently deter mine the distribution o f rest-mass and energy-momentu m injected, 
which is left for future work. In this paper, the rough results of Popham et al. ( 1999); MacFadven & Wooslev ( 1999) are used to 
approximate the neutrino annihilation and pair creation. They determine the energy density creation rate in pairs as measured by 
an observer at infinity (<*,„,). A monoenergetic, monomomentum injection of rest-mass and energy-momentum is assumed. That 
is, for a Lorentz invariant particle distribution function of 

/= ., (A4) 

dx l dx 2 dx 3 du 1 diP-du? 

the fluid equations are derived from the Boltzmann equation 

df dn 
dr 



(A5) 



for a particle creation density rate dn^j/dr in the comoving frame of the existing fluid. There is no collisional term in the ideal 
case. After mass-weighting and taking 4-velocity moments of the Boltzmann equation one has S p = (po,inj u ^,j);tM> where po./nj is 
the injected rest-mass density and uf n j is the 4-velocity of the injected particles. Also, Sj = T/%" . Here, 



dn h 



mj 



m e dn, 



e e + ,mj 



1 dpo,e-e+jnj 



dr 



nit, dr 



m h 



dr 



(A6) 



18 



McKinney, J.C. 



so the fluid is still treated as a single component with a single mass weight mi,. As discussed in iMcKinnevI (l2005bl) . this is a 
good approximation since the late-time Poynting-dominated jet region is dominated by lepton mass while the remaining region 
is dominated by baryon mass. The simultaneous advected flux of injected energy is neglected, so then St, = J—ge-mj is valid in 
Boyer-Lindquist, Kerr-Schild, or modified Kerr-Schild coordinates. This two-step approach (similar to operator splitting), where 
particles are injected and then advected with the normal fluid, is a good approximation for typical numerical integrations that use 
a multi-step timestep approach. This is also a good approximation because the injection region has a small 4-velocity beyond the 
stagnation surface. Also, the error in the injection approximations being made is larger than the error in neglecting the advection 
term. Thus 

S p = d/dt(y^gp^ inj u' inj ), (A7) 

and 

St, = d/dt y=g {qojnj U ' inj ^ i+ P'M) ■ < Ag ) 
where qo,inj = po,inj+Uinj+Pinj- This represents the source of energy-momentum, such as photon and neutrino losses emitted and 
absorbed and the pairs created by annihilation of photons or neutrinos. This accounts for, in the guise of the MHD approximation, 
"non-local" transport of energy and momentum. The injected particles are a relativistic gas with 7 = 4/3 such that p m j = (7- 1 )«;«/. 

Two approximate injection methods are used to treat the injected momentum. One assumes that the momentum injected is 
mostly aligned with the axis with uf n j = and either vf n j = or vf n j = il[R = r stag ] for some typical origin of the neutrinos r stag . 
For the other, the particles are injected in the comoving frame of the existing fluid (i.e. u^ n - = u^). It turns out that these different 
approaches lead to qualitatively similar results for the solution beyond the region where injection is important (r > 6ri in the 
jet). The total energy density injected (e',„ ; as defined by an observer at infinity) is partitioned between rest-mass, enthalpy, and 
momentum energy. The fraction of rest-mass injected is defined as f p and so 

d/dt(y^gpo M ju' inj ) = f p ^/=ge in j. (A9) 
The fraction of internal energy related injected energy is defined as //,, and so 

9/ dt{^g((u in j + p in j)u' inj U t nJ + Pin j)) = f h yf^iin j. (A 1 0) 

The fraction of rest-mass plus momentum energy injected is defined as f p +f m , and so 

d/dt(y/=gp ,i nj u , inju' l nj ) = (fp+fmlV^geinj- (All) 

Here f n +h+ f„= 1 

McKinnev (2005b) defines e\ n j, which along with the injection fractions and u^u p = -1 completely define the injection process 
by allowing one to solve for p m j, Ui„j, u\ n p and u'" J for a fixed time interval dt. Radiative cooling effects are important for studying 
hundreds of dynamical times, but otherwise can be neglected if the initial model is consistent with the time-averaged radiative 
properties. Future work will refine the treatment of the pair creation process, considering the creation effects in the locally flat 
comoving frame, once a self-consistent model is available. 

A.2. Pair Plasma Annihilation 

The electron-positron pair plasma that forms may annihilate itself into a fireball if the pair annihilation rate is faster than 
the typical rate of the jet (c 3 /GM) near the black hole. Also, if the pair annihilation timescale is shorter than the dynamical 
ti me, then pa ir annihilation would give a collisional term in the Boltzmann equation. From the pair annihilation rate given 
m McKinney (2005b), one finds that t pa GM/c 3 for AGN and marginally so for x-ray binaries. Thus, pairs mostly do not 
annihilate, and so formally the pair plasma that forms in the low-density funnel region is collisionless so that the Boltzmann 
equation should be solved directly. Plasma instabilities and relativistic collisionless shocks are implicitly assumed to keep the 
pairs in thermal equilibrium so the fluid approximation remains mostly valid, as is a good approximation for the solar wind (see, 
e.p. lFeldman & Marschll997tlUsmanov et al.l2000|). This same ap proximation has to be invoked for the thick disk state in AGN 
and x-ray binaries, such as for the ADAF model dMcKinne\l2004l) . For regions that pair produce slower than the jet dynamical 
time, each pair-filled fluid element has a temperature distribution that gives an equation of state with P = po,e-<&kbT e /m e rather 
than P = (ll/12)ar 4 , where a is the radiation constant. So most of the particles have a Lorentz factor of T e ~ u/(po e - e +c 2 ) 
and little of the internal energy injected is put into radiation. This also allow the use of a single-component approximation. A 
self-consistent Boltzmann transport solution is left for future work. 

On the contrary for GRB systems, due to the relatively high density of pairs, the time scale for pair annihilation is t pa <C GM/c 3 
along the entire length of the jet. Thus a pair fireball forms and the appropriate equation of state is that of an electron-positron- 
radiation fireball. Thus, formally the pair fireball rest-mass density is not independent of the pair fireball internal energy density. 
However, because the pairs are well-coupled to the radiation until a much larger radius of r ~ 10 8 - 10 10 tl, the radiation provides 
an inertial drag on the remaining pair plasma. That is, the relativistic fluid energy-momentum equation is still accurate. So the 
effective rest-mass density is ~ po + u (u the total internal energy of the firebal l), and so the e ffective rest-mass is independent of 
the cooling of the fireball until the fireball is optically thin (see, e.g., Meszaros & Rees 1997). 

For GRB systems, the mass conservation equation is reasonably accurate. Even though the electron-positron pairs annihilate, 
th e rest-mass of pairs injected is approximately that of the pairs that are injected due to Fick-diffusion of neutrons (see appendix A 
of McKinnev 2005b). The annihilation energy from electron-positron pairs contributes a negligible additional amount of internal 
energy, so can be neglected, especially compared to the Poynting energy flux that emerges from the black hole. Thus, the rest- 
mass can always be assumed to be due to baryons rather than the electron-positron pairs. This also suggests that the neutrino 
annihilation is a negligible effect if the BZ power is larger than the neutrino annihilation power. 

In summary, the rest-mass evolution discussed in section|3]is accurate for GRB, AGN, and marginally so for x-ray binaries. 
This is despite the lack of Boltzmann transport for the collisional system, or a collisional term due to pair annihilation. 
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A.3. Electromagnetic Terms 
In terms of F^ u , the Faraday (or electromagnetic field) tensor, 



TZ=F^F%-^F a f 3 F a0 , (A12) 



where a factor of v4tt is absorbed into the definition of . The induction equation is given by the space components of 
F ; „ = 0, where F is the dual of the Faraday, and the time component gives the no-monopoles constraint. The other Maxwell 
equations, = F^ v define the current density 7 M . The comoving electric field is defined as 

e v = = l -e^ kx u u *F xk = rjf, (A13) 

where 77 corresponds to a scalar resistivity for a comoving current density = J V P U ^, where P v,i = + u v u^ projects any 
4-vector into the comoving frame (i.e. P vy, u^ = 0). The classical ideal MHD approximation that r\ = e^ = is assumed. The 
comoving magnetic field is defined as 

V = u*F^ = U^ kX u u *F xk , (A14) 



and so the stress-energy tensor can be written as 

b 2 
2 

and 



S = -( bV + H-W, (A15) 



and so 



V* ' = W-iV (A16) 



(A17) 



since F = ^e^ VKX F K \. Here e is the Levi-Civita tensor. Following the notation of MTW, e^ xs = -^^[fivXS], where [pvAS] 

is the completely antisymmetric symbol and = 0, 1, or -1. Notice that e v u v = b v u v = 0, so they each have only 3 independent 
components and are space-like 4- vectors. 

With B' = F and E' = F", the no-monopoles constraint becomes 

(^=gB% = 0, (A18) 

and the magnetic induction equation becomes 

(V=gB% = -(^g{b i u i -bij)), ] . (A19) 

The ideal MHD approximation assumes that e M = 0, and so the invariant e^b^ = 0. Since the Lorentz acceleration o n a pa rticle is 
fl = qe^, then this implies that the Lorentz force vanishes on a particle in the ideal MHD approximation. Equation lA14l implies 
b' = B'lij and V = (B' + u'b') /«', so the magnetic induction equation becomes 

( = -(^(fi V - W) j 

=-(V^e y *e*)),y, (A20) 

where v ! = u'/u', e, = —ei^B k = -v x B is the E MF, and t l,k is the spatial permutation tensor. A more complete account of the 
relativistic MHD equations can be found in Anile] ( I19891) . 

B. CHARACTERISTIC (AND OTHER) SURFACES 

The ideal MHD dispersion relation is given, e.g., in lGammie et alJ J200 3a) (there is a sign typo there), and summarized here. 
In the comoving frame, the dispersion relation is 

D(k») = 

Lu(LU 2 -(k-V A ) 2 ) X (Bl) 

(lu 4 -lu 2 (K 2 c 2 ns + c 2 (k-v A ) 2 /c 2 )+K 2 c 2 (k-v A ) 2 ) , 

where c 2 ns = (\\ + c 2 (l - \ 2 A /c 2 )) is the magnetosonic speed, c 2 = (d(p+u)/dp)~ x is the relativistic sound speed, Va = B/y£ is the 
relativistic Alfven velocity, £ = b 2 + w, and w = p+u + p. Here c is the (temporarily reintroduced) speed of light. The invariant 
scalars defining the comoving dispersion relation are uj = k^u^, K 2 = K^K* = k^+uj 2 , where = P^k" = k^ + wm^ is the 

projected wave vector normal to the fluid 4-velocity, v A = b^b^/S, and (k • va) = k^b^/VT. The terms in the dispersion relation 
correspond to, respectively from left to right, the zero frequency entropy mode, the left and right going Alfven modes, and the 
left and right going fast and slow modes. The eighth mode is eliminated by the no-monopoles constraint. 

The dispersion relation gives the ingoing and outgoing slow, Al fven, and fast surface s. Energy can be extracted from the 
black hole if and only if the Alfven point lies inside the ergosphere (Takahashi et al. 1990). Opti mal acceleratio n of the flow by 
conversion of Poynting flux to kinetic energy flux occurs beyond the outer fast surface dBeeelman Other surfaces 
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include: the horizon at r H = 1 + yj 1 - j 2 ; the ergosphere at r = 1 + ^/l-(jcos#) 2 ; the coordinate basis light surface in where 
^/g^p = c/Qf, where asymptotically ^/g^p = rsin(fJ) is the Minkowski cylindrical radius ; the surface in Boyer-Lindquist coor- 
dinates where the toroidal field equals the poloidal field where = B r . Finally, there is a stagnation surface where the poloidal 
velocity u p = 0. In a field confined jet where no matter can cross field lines into the jet, and if the jet has inflow near the black 
hole and outflow far from the black hole, then this necessarily marks at least one location where rest-mass must be created either 
by charge starving the magnetosphere till the Goldreich-Julian charge density is reached, or pair production rates sustains the 
rest-mass density. 
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